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Abstract 

D I We present updated leading-order, next-to-leading order and next-to-next-to-leading 

j-^' order parton distribution functions ("MSTW 2008") determined from global analysis of 

hard-scattering data within the standard framework of leading-twist fixed-order collinear 
I factorisation in the MS scheme. These parton distributions supersede the previously avail- 

^ ■ able "MRST" sets and should be used for the first LHC data-taking and for the associated 

. theoretical calculations. New data sets fitted include CCFR/NuTeV dimuon cross sections, 

which constrain the strange quark and antiquark distributions, and Tevatron Run II data 
on inclusive jet production, the lepton charge asymmetry from W decays and the Z rapid- 
ity distribution. Uncertainties are propagated from the experimental errors on the fitted 
I data points using a new dynamic procedure for each eigenvector of the covariance matrix. 

' We discuss the major changes compared to previous MRST fits, briefly compare to parton 

distributions obtained by other fitting groups, and give predictions for the W and Z total 
^ cross sections at the Tevatron and LHC. 
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1 Introduction 



In deep-inelastic scattering (DIS), and "hard" proton-proton (or proton-antiproton) high-energy 
coUisions, the scattering proceeds via the partonic constituents of the hadron. To predict the 
rates of the various processes a set of universal parton distribution functions (PDFs) is required. 
These distributions are best determined by global fits to all the available DIS and related hard- 
scattering data. The fits can be performed at leading-order (LO), next-to- leading order (NLO) or 
at next-to-next-to- leading order (NNLO) in the strong coupling as- Over the last couple of years 
there has been a considerable improvement in the precision, and in the kinematic range, of the 
experimental measurements for many of these processes, as well as new types of data becoming 
available. In addition, there have been valuable theoretical developments, which increase the 
reliability of the global analyses. It is therefore timely, particularly in view of the forthcoming 
experiments at the Large Hadron Collider (LHC) at CERN, to perform new global analyses 
which incorporate all of these improvements. 

The year 2008 marked the twentieth anniversary of the publication of the first MRS analysis 
of parton distributions [H, which was, indeed, the first NLO global analysis. Following this 
initial MRS analysis there have been numerous updates over the years, necessitated by both the 
regular appearance of new data sets and by theoretical developments |2h21[|. In the modern era, 
the MRST98 (NLO) sets 12| were the first of our sets to take full advantage of a large amount 
of new HERA structure function data, and the first to incorporate heavy quarks in a consistent 
and rigorous way in the default PDF set (using a general-mass variable flavour number scheme 
approach; see Section Hj). The uncertainty in the PDFs was explored by producing a modest 
number of additional sets with different gluon distributions and different values of the strong 
coupling. The following year, the NLO set was updated (MRST99), a LO set was produced, 
and new sets corresponding to a different treatment of higher-twist contributions and the DIS 
scheme were presented |131]. The year 2001 saw a major update [14 1, w ith a NNLO (MRST 2001) 
set, based on an approximation of the NNLO splitting functions 2a, l23|], produced for the first 



time [15|, alongside new NLO and LO sets. The grid interpolation was also improved to allow 
for faster and more accurate access to the PDFs in the public interface code. The following 
year the Hessian approach (see Section [H]) was used to produce a "parton distributions with 
errors" package (MRST 2001 E) comprising a central NLO set and 30 extremum sets (l6l| . The 
central NLO and NNLO sets were updated slightly in the same year (MRST 2002) [l6|, using 
an improved approximation to the NNLO splitting functions 2J|. In 2003, fits were performed 
in which the x and range of DIS structure function data was restricted to ensure stability 
with respect to cuts on the data, and corresponding NLO and NNLO "conservative" variants 
of the MRST 2002 sets were derived (MRST 2003 C) [l7|. The next major milestone was in 
2004, with a substantial update of the NLO and NNLO sets (MRST 2004) [ISj], the latter using 
the full NNLO splitting functions 25|, l26| for the first time and both incorporating a "physical" 
parameterisation of the gluon distribution in order to better describe the high-ET Tevatron jet 
data. A NLO set incorporating 0{a) QED corrections in the DGLAP evolution equations was 
also produced for the first time (MRST 2004 QED) [l9| , together with fixed-fiavour-number LO 
and NLO variants [i^]. Finally, in 2006 a NNLO set "with errors" was produced for the first 
time (MRST 2006 NNLO) j2l[, using a new general-mass variable flavour number scheme and 
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MSTW 2008 NLO PDFs (68% C.L.) 
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Figure 1: MSTW 2008 NLO PDFs at = 10 GeV^ and = 10^ GeV^ 



with broader grid coverage in x and than in previous sets. 

In this paper we present the new MSTW 2008 PDFs at LO, NLO and NNLO. These sets are 



a major update to the currently available MRST 2001 LO [15|, MRST 2004 NLO [18( and MRST 



2006 NNLO [21] PDFs. The "end products" of t he p resent paper are grids and interpolation 
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An example is given in Fig. [T|, which 
10 GeV^ and = 10^ GeV^, including the associated 



code for the PDFs, which can be found at Ref. 
shows the NLO PDFs at scales of 
one-sigma (68%) confidence level (C.L.) uncertainty bands. 

The contents of this paper are as follows. The new experimental information is summarised in 
Section [2J An overview of the theoretical framework is presented in Section [3] and the treatment 
of heavy flavours is explained in Section HI In Section Owe present the results of the global flts and 
in Section [6] we explain the improvements made in the error propagation of the experimental data 
to the PDF uncertainties, and their consequences. Then we present a more detailed discussion of 
the description of different data sets included in the global flt: inclusive DIS structure functions 
(Section [7]), dimuon cross sections from neutrino-nucleon scattering (Section [H]), heavy flavour 
DIS structure functions (Section [H]), low-energy Drell-Yan production (Section [TD]), W and Z 
production at the Tevatron (Section [11]), and inclusive jet production at the Tevatron and 
at HERA (Section [T2l) . In Section [13] we discuss the low-x gluon and the description of the 
longitudinal structure function, in Section [H] we compare our PDFs with other recent sets, 
and in Section [TJ] we present predictions for W and Z total cross sections at the Tevatron and 
LHC. Finally, we conclude in Section [TH] Throughout the text we will highlight the numerous 
refinements and improvements made to the previous MRST analyses. 
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2 Survey of experimental developments 



Since the most recent MRST analyses [15|, Il8|, l21[ a large number of new data sets suitable for 
inclusion in the global fit have become available, or are included for the first time. Some of 
these are entirely new types of data, while others supersede existing sets, either improving the 
precision, extending the kinematic range, or both. Here, we list the new data that we include in 
the global fit, together with an indication of the parton distributions that they mainly constrain. 



Compared to the analysis in Ref. there is no new large ep — > eX deep-inelastic F2 
structure function data set. However, note that we now fit to the measured reduced cross 
section values 

a{x, Q') = F,{x, Q') - ^—^—^^Fl{x, Q'), (1) 
where y = Q'^Jixs) and ^/s is the ep centre-of-mass energy. In fact, this was already done 



in Refs. |2ll . |28| . Since y -C 1 in most of the kinematic range, a{x,Q^) is effectively the 



same as F2{x, Q^). However, at HERA, for the lowest x- values at a given Q^, the value of 



30|. 



y can become as large as 0.7-0.8, and the effect of Fl{x,Q^) becomes apparent [29 
We now include a small, but important, number of data points from HI in Ref. 
The effect of Q^) is seen in the data as a flattening of the growth of a{x,Q^) as 

X decreases to very small values (for fixed Q^), leading eventually to a turnover. Hence, 
for precise analysis of the HERA structure function data it is particularly important to 
fit any theoretical prediction to the measured a{x,Q'^), rather than to model-dependent 
extracted values of F2{x, Q^). In previous analyses we have used the next best approach to 
fitting to a{x, Q^) directly by making the correction using our own prediction for Fl{x, Q^), 
maintaining self-consistency. We now also include data on Fl{x, Q^) at high x from fixed- 



target experiments j32H34| . which provide a weak constraint on the gluon for x ~ 0.1. 



(ii) The structure functions F2 and xF^ for the charged-current deep-inelastic processes e p 



UeX and e p UeX have been measured at HERA [35|, l36l]. At large x these processes are 
dominated by the u and d valence quark distributions respectively, although the event rate 
is low in this domain. In principle, these data are valuable as they constrain exactly the 
same partonic combinations as the data which exist for the neutrino processes u^N fi'^X 
and Ufj^N fJ''~X, but without the problems of having to allow for the effects of the heavy 
nuclear target, A^, but the precision of published data is still relatively low. We include 
the data on e^p UeX since these are more precise, and constrain the less well-known 
down valence distribution. 



(iii) We include new neutrino structure function data from NuTeV (with an iron target) [3j 
and CHORUS (with a lead target) jssj, though consistency requires the use of only data 
with X < 0.5, as discussed later in Section These provide information on the flavour 
decomposition of quarks for 0.01 < x < 0.5. 

(iv) Data also from the neutrino DIS experiments NuTeV and CCFR on dimuon production 



give a direct 39| constraint on the strange quark content of the proton for 0.01<x<0.2, 
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VI 



IX 



replacing the previous direct, but mo del- dependent, constraint on the strange distribu- 
tion ^o|. Furthermore, the separation into neutrino and antineutrino cross sections allows 
a separation into quark and antiquark contributions and therefore a determination of the 
strangeness asymmetry. 



There are now improved measurements of the structure functions, Fg'^ and 4lN47 



for heavy quark production at HERA. We include all of the published data on F^'^ in the 
analysis. Since these HERA data are driven by the small-x gluon distribution, via the 
g HH transition, they probe the gluon for x ~ 0.001. However, they also provide useful 
information on the mass of the charm quark and are a test of our procedure for including 
heavy flavours. Since the data on are far less numerous and of lower precision, they 
do not provide any further constraints on the gluon and so we simply compare to these. 



There are Tevatron Run II data on the lepton charge asymmetry from W decays [48H50 . 
These are more precise than the Run I data |51| and provide constraints on the down quark 
distribution for x > 0.05, and also to a smaller extent on the up and down sea quarks for 



X 



0.1. 



(vii) For the first time there are data from CDF 52| and DO 53| on the Z (— s> l^l ) rapidity 



distribution which constrain the quarks for x > 0.05 with a different weighting than DIS 
data and at a significantly higher scale. 

viii) There now exist data for inclusive jet production from Run II at the Tevatron from 



CDF |54j . l55l | and D0 [56|]. These data exist in different rapidity bins, are more precise 
than the previous Run I data and go out to larger jet pt values. They provide constraints 
on the gluon (and quark) distributions in the domain 0.01 < a; < 0.5. 

We include data from HERA on inclusive jet production in DIS 57-5^. (We do not include 
the photoproduction data because of potential sensitivity to photon PDFs.) These data 
constrain the gluon for 0.01 < x < 0.1. 

The total data that we use in the new global analyses are listed in Table [2] of Section [5], 
together with the individual values for each data set for the LO, NLO and NNLO fits. A 
rough indication of the particular parton distributions that the various data constrain is also 
given in Table 1. 
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Process Subprocess Partons x range 

£^ {p, n} ^ X 7*g — >• q q, q, g x > 0.01 

e^n/p^e^X Yd/u^d/u d/u x > 0.01 

pp fx+fj,- X uM.Jd ^7* q_ 0.015 < x < 0.35 

pn / pp ^ fi~^ fi^ X {ud)/{uu)^j* d/u 0.015 < x < 0.35 

u{iy)N ^ ^^-{^l+)X W*q^q' q,q 0.01 < x < 0.5 

uN^fi-fj.+ X W*s^c s 0.01 < a; < 0.2 

VN ^ ji+ji-X W*s^c I 0-01 < a: < 0.2 

e±p^e±X 7*9^9 0.0001 < x < 0.1 

e+p^PX iy+ {d, s} ^ {m, c} ci,s a; > 0.01 

e^p ccX 7*c c, 7*5f ^ cc c, 51 0.0001 < .x < 0.01 

e^p ^ jet + X qq g o.oi < x < o.i 

pp ^ iet + X gg, qg, qq 2j g,q _ 0.01 < a; < 0.5 

pp {W^ e^u) X ud^W,ud^W u,d,u,d x> 0.05 

pp^ {Z ^ e+e-)X uu,dd^Z d x>0.05 



Table 1: The main processes included in the current global PDF analysis ordered in three groups: 

fixed-target experiments, HERA and the Tevatron. For each process we give an indication of their 
dominant partonic subprocesscs, the primary partons which are probed and the approximate 
range of x constrained by the data. 
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3 Overview of theoretical framework 



In this section we first give a brief overview of the standard theoretical formahsm used, and then 
present a summary of the theoretical improvements and changes in methodology in the global 
analysis. A more detailed discussion of the various items is given later in separate sections. 

We work within the standard framework of leading-twist fixed-order coUinear factorisation 
in the MS scheme, where structure functions in DIS, Fi{x,Q'^), can be written as a convolution 
of coefficient functions, Ci^a, with PDFs of flavour a in a hadron of type A, fa/A{x, Q^), i.e. 

Fi{x, = ^ ® fa/A{x, Q'). (2) 

a=q,g 

Similarly, in hadron-hadron collisions, hadronic cross sections can be written as pro cess- dependent 
partonic cross sections convoluted with the same universal PDFs, i.e. 

O-AB = ^ O-afe ® fa/Aixi, Q"^) ® fb/B{x2, Q^)- (3) 
a,b=q,g 

The scale dependence of the PDFs is given by the DGLAP evolution equation in terms of 
the calculable splitting functions, Paa', i-e. 

m^-^P'-^f-''^^ (4) 

a'=q,g 



The DIS coefficient functions, Ci^a, the partonic cross sections, aab, and the splitting functions, 
Paa', can each be expanded as perturbative series in the running strong coupling, asiQ"^)- The 
strong coupling satisfies the renormalisation group equation, which up to NNLO reads 

^\ = -pJ^Y^gJ^Y-^Jm'-.... (5) 



dlnQ2 V47r/ KAttJ KAtc J KAtc J 

The input for the evolution equations, (jl]) and fa/A{x,Ql) and asiQo)^ ^ reference 
input scale, taken to he Ql = 1 GeV^, must be determined from a global analysis of data. In 
the present study we use a slightly extended form, compared to previous MRST fits, of the 
parameterisation of the parton distributions at the input scale Ql = 1 GeV^: 



XUy (X, Qq) 


= A^x'^-{l- 


- + e„ a/x + -fux), 


(6) 


xd^{x^ Qq) 


= AdX^-^{l- 


-xY'^il + edV^ + ^dx), 


(7) 


xS{x,Ql) 


= Asx'^{l^ 


-xy^{l + esV^ + lsx), 


(8) 


x/^[x,Ql) 


= Aax''^(1 


-x)''s+2(1+7ax + 5ax'), 


(9) 


xg{x,Ql) 


= Agx^^{l - 


- X^'il + eg ^/x + -fg X) + Ag, X ^ ( 1 " X) V , 


(10) 


x{s + s){x,Ql) 


= A+ x^^ (1 


-xy+{l + esV^ + ^sx), 


(11) 


x{s - s){x,Ql) 


= A_ x^-{l 


- x)''-(l - x/xo). 


(12) 
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where A = d — u, = q — q, and where the hght quark sea contribution is defined as 

S = 2{u + d) + s + s. (13) 
The input PDFs hsted in Eqs. fl6l)- f|T2|) are subject to three constraints from number sum rules: 

/ dx Uy{x,Ql) = 2, / dx dy{x,Ql) = 1, / dx Sy{x,Ql) = 0, (14) 
Jo Jo Jo 

together with the momentum sum rule: 

/ dx X [uy{x,Ql) + d^{x,Ql) + S{x,Ql) + g{x,Ql)] = 1. (15) 
Jo 

We use these four constraints to determine Ag, A^, A^ and Xq in terms of the other parameters. 
There are therefore potentially 34 — 4 = 30 free PDF parameters in the fit, including as- The 
values of the parameters obtained in the LO, NLO and NNLO fits are given in Table H] in 
Section [5] below. (In practice, we fix 5- = 0.2 in (|T2l) due to extreme correlation with A^ and 
?7_.) For the LO fit, where there is no tendency for the input gluon distribution to go negative 
at small x, the second term of the parameterisation ffTOl) is omitted. 

The major changes to the analysis framework, compared to previous MRST analyses, are 
listed below. 

(i) We produce and present PDF sets at LO, NLO and NNLO in as- The last of these uses 
the splitting functions calculated in Refs. (25|^l26^ and the (massless) coefficient functions 
for structure functions calculated in Refs. j60h65.] together with the massive coefficient 
functions described in Section HI We do not include here a set of modified LO PDFs for 



use in LO Monte Carlo generators of the form described in Refs. [66|, |67l]. This will be the 
topic of a further study. We produce eigenvector sets describing the uncertainty due to 
experimental errors in all cases. This is the first time this has been done at LO. We have 
significantly modified the method of determining the size of the uncertainties, no longer 
using a simple Ax^ = 50 above the global minimum to determine the 90% confidence level 
uncertainty for the PDFs. We now have a dynamic determination of the tolerance for each 
eigenvector direction, which gives qualitatively similar results to previously. That is, we 
now have no a priori fixed value for the tolerance. This method is described in detail in 
Section [6l We also include, for the first time, the uncertainties in the PDFs due to the 
uncertainties in the normalisations of the individual data sets included in the global fit. We 
present the grids for the eigenvector PDFs as both one-sigma (68%) and 90% confidence 
level uncertainties, rather than just the latter. 



I.e. we 



(ii) We now use the simpler and more conventional definition of defined in Ref. |68l ], 

solve ([5]) truncated at the appropriate order, starting from an input value of as{Ql)- 
This is particularly important at NNLO where the coupling becomes discontinuous at the 
heavy fiavour matching points. Further discussion of the difference between the MRST 
and MSTW definitions of as will be given elsewhere. The value of as{Ql) is now one 
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of the fit parameters and replaces the Aqcd parameter used in MRST fits. Since there 
is more than one definition of Aqcd in common use, replacing it by as{Qo) reduces the 
potential scope for misuse. 

(iii) The new definition of the coupling is then identical to that used in the PEGASUS jisl 



(which uses moment space evolution) and hoppet [69| (which uses x space evolution) 
programs, therefore we can now test our evolved PDFs by benchmarking against those 
two programs, which were shown to display excellent agreement with each other 70. 7ll . 



We use our best-fit input parameterisations at Ql = 1 GeV^ in each case, which are more 
complicated than those in the toy PDFs used for the benchmarks in Refs. LDJ, and 
perform trial evolutions at LO, NLO and NNLO up to = 10^ GeV^, finding agreement 
with both PEGASUS and hoppet to an accuracy of (9(0.1%) or less in most regions, with 
discrepancies a little larger only in regions where the PDFs are very small or at extremely 
low X where extrapolations are used. (We see larger discrepancies with PEGASUS at a; ~ 0.9, 
where the default choice of the Mellin-inversion contour is inaccurate when using our input 
parameterisation.) We intend to return to more comprehensive comparisons in the future, 
but for the moment are satisfied that any discrepancies between evolution codes are orders 
of magnitude lower than the uncertainties on the PDFs. 

(iv) For structure function calculations, both neutral and charged current, we use an improved 
general-mass variable flavour number scheme (GM-VFNS) both at NLO and more partic- 



ularly at NNLO. This is based on the procedure defined in Ref. [72], though there are some 
slight modifications described in detail in Section |H The new procedure is much easier to 
generalise to higher orders. At NLO there is simply a change in the details, but at NNLO 
there is a major correction to the transition across the heavy flavour matching points. 



first implemented in Ref. [21|, which demonstrates that NNLO sets before 2006 should 
be considered out-of-date. No other available NNLO sets which include heavy flavours 
currently treat the transition across the matching points in such a complete manner. A 
full GM-VFNS is not currently available for Drell-Yan production of virtual photons, or 
W and Z bosons, or for inclusive jet production. For the low-mass Drell-Yan data in our 
flt, heavy flavours contribute <^ 1% of the total, so the inaccuracy invoked by the ap- 
proximation of using the zero-mass scheme is negligible. All other processes are at scales 
such that the charm mass rric is effectively very small, and the zero-mass scheme is a very 
good approximation. The approximation should still be reliable for processes induced by 
bottom quarks, and in this case the relative contribution is small. 

Vector boson production data can now be described in a fully exclusive way at NNLO. We 



use the FEWZ code [73|, and compare with the results obtained using the resbos code [74 
which includes NLO-I-NNLL p^^-resummation effects. This allows, for example, a detailed 
description of the W asymmetry data accounting for the W width (Tiy) and lepton decay 
effects. 



vi) We now implement fastnlo [75||, based on nlojet-| — \- [76j, |77j, which allows the inclusion 



of the NLO hard cross section corrections to both the Tevatron and HERA jet data in the 
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fitting routine. This improvement replaces the i^-factors and pseudo-gluon data that were 
previously necessary to speed up the fitting procedure. 



(vii) We provide [27[ a new form of the grids that list the parton distributions over an extended 
range and which are more dense in x. The extended public grids allows an improved 
description around the heavy quark thresholds, where at NNLO discontinuities in parton 
distributions appear. (We also make improvements to the internal grids used in our evolu- 
tion code, so that the heavy quark thresholds lie exactly on the grid points, which was not 
the case in the MRST analyses.) Moreover, the grids now contain more parton flavours: 
s 7^ s at all orders due to the (mild) evidence from the dimuon data, and additionally 
c 7^ c and 6 7^ 6 at NNLO, due to automatic generation of a small C(a|) quark-antiquark 
difference during the NNLO evolution, so at this order the H distribution is (very) slightly 



different from H 781]. (This small asymmetry at NNLO was omitted in previous MRST 
analyses.) 



, we 



(viii) Two features of the parameterisation are worth emphasising. First, as since Ref. [14 

allow the gluon to have a very general form at low x, where it is by far the dominant 
parton distribution. We will discuss the consequences of this in Section 16.51 Second, we 
parameterise s ± s as functions of x, rather than assuming, as previously, 

s = s = -{u + d) (16) 

at the input scale, where k was a constant fixed by the neutrino dimuon production (/t ~ 
0.4-0.5); see Section [8l Implicit in our parameterisation is the assumption that the strange 
sea will have approximately the same shape as the up and down sea at small x. This extra 
parametric freedom of the gluon and strange quark distributions means that all partons 
are less constrained, particularly at low x. This, in turn, leads to a more realistic estimate 
of the uncertainties on all parton distributions. 

(ix) We fix the heavy quark masses at rric = 1.40 GeV and = 4.75 GeV, changed from the 
MRST default values of = 1.43 GeV and mf, = 4.30 GeV. Our value of = 4.75 GeV 
is close to the calculated MS mass transformed to the pole mass value, using the three-loop 
relation between the pole and MS masses, of = 4.800 GeV 79, 80], while our value of 
rric = 1.40 GeV is smaller than the calculated pole mass value of rric = 1.666 GeV 79, 8^ . 



If allowed to go free in the global fits, the best-fit values are nic = 1.39 GeV at NLO and 
rric = 1.27 GeV at NNLO. In a subsequent paper we will present a more detailed discussion 
of the sensitivity of different data sets to the charm quark mass and present the best-fit 
values including a determination of the uncertainty. We allow a maximum of five flavours 
in the evolution and do not include top quarks. 
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4 Treatment of heavy flavours 



The correct treatment of heavy flavours in an analysis of parton distributions is essential for 
precision measurements at hadron colliders. For example, the cross section for W production 
at the LHC depends crucially on precise knowledge of the charm quark distribution. Moreover, 
a correct method of fitting the heavy flavour contribution to structure functions is important 
because of the knock-on effect it can have on all parton distributions. However, it has become 
clear in recent years that it is a delicate issue to obtain a proper treatment of heavy flavours. 
There are various choices that can be made, and also many ways in which subtle mistakes can 
occur. Both the choices and the mistakes can lead to changes in parton distributions which may 
be similar to, or even greater than, the quoted uncertainties — though the mistakes usually lead 
to the more dramatic changes. Hence, we will here provide a full description of our procedure, 
along with a comparison to alternatives and some illustrations of pitfalls which must be avoided. 

First we describe the two distinct regimes for describing heavy quarks where the pictures are 
relatively simple. These are the so-called fixed flavour number scheme (FFNS) and zero-mass 
variable flavour number scheme (ZM-VFNS); see Fig. [21 



Figure 2: Diagrams contributing at leading-order to the heavy flavour structure function, , in 
(a) the fixed flavour number scheme (FFNS) valid for < m|^, and (b) the zero-mass variable 
flavour number scheme (ZM-VFNS) valid for ^ m'jj. 



4.1 Fixed flavour number scheme 

First, there is the region where the hard scale of the process is similar to, or smaller than, the 
quark mas^, i.e. < mj^. In this case it is most natural to describe the massive quarks as 
final-state particles, and not as partons within the proton. This requirement defines the FFNS, 
where only light quarks are partons, and the number of fiavours is fixed. We label the number of 
quark flavours appearing as parton distributions hj Uf. Note, however, that there are multiple 
instances of FFNSs with different numbers of active flavours Uf. The number rif is normally 
equal to 3, where up, down and strange are the light quarks, but we can treat charm as a light 

"'^Throughout this section we use H to denote a heavy quark: H ^ c, b or t. 
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quark while bottom remains heavy at high scales, i.e. Uf = 4, and it may also be equal to 5 if 
only top quarks are treated as heavy. In each example of a FFNS the structure functions are 
given bj0 

F.(x,g^) = $:<r-^(gVO ® fk'm- (17) 



This approach contains all the mj^-dependent contributions, and as it is conceptually simple it is 
frequently used in analyses of structure functions. Even in this case one must be careful to be self- 
consistent in defining all quantities, i.e. parton distributions, coefficient functions and coupling 
constant in the same renormalisation schemes (which is often not done). The mistake mad e by 
not doing so can lead to errors in the gluon distribution of order the size of the uncertainty 20| . 

Despite its conceptual simplicity, the FFNS has some potential problems. It does not 
sum W (Q'^ / TTijj) {I < m) terms in the perturbative expansion. Thus the accuracy of the 
fixed-order expansion becomes increasingly uncertain as increases above rn^j. For example, 
0^5 (M^) \n{M^/ml) is approximately equal to unity, and so there is no guarantee an expansion 
in this variable will converge. As well as this problem of principle, there are additional practical 
issues. Since calculations including full mass dependence are complicated, there are only a few 
cross sections known even to NLO in as within this framework, so the resulting parton distri- 
butions are not universally useful. In addition, even for neutral-current structure functions, the 



FFNS coefficient functions are known only up to NLO [8l|, |82|, and are not calculated at NNLO 
— that is, the a| coefficient!!, C2^^'^'^^\ for F2 is unknown, so one cannot determine parton 
distributions at NNLO in this scheme. 



4.2 Zero-mass variable flavour number scheme 

All of the problems of the FFNS are solved in the so-called zero-mass variable flavour number 
scheme (ZM-VFNS). Here, the heavy quarks evolve according to the splitting functions for 
massless quarks and the resummation of the large logarithms in Q'^/m'jj is achieved by the 
introduction of heavy-flavour parton distributions and the solution of the evolution equations. 
It assumes that at high scales, ^ m^, the massive quarks behave like massless partons, and 
the coefficient functions are simply those in the massless limit, e.g. for structure functions 

F.(x,Q^) = 5^C™'"^®/;^(Q^), (18) 

k 

where Uf — 3 is the number of active heavy quarks, with masses above some transition point for 
turning on the heavy flavour distribution, typically at a scale similar to m|^. This is technically 
simpler than the FFNS, and many more cross sections are known in this scheme. The nomen- 
clature of "zero-mass" is a little misleading because some mass dependence is included in the 

^To simplify the discussion, we will use the convention that the factorisation scale and renormalisation scale 
are both equal to Q^. It is also the choice we make in the analysis. Alternative choices are possible, and cause 
no problems in principle, but can add considerable technical complications. 

■^We add a subscript H to distinguish the g ^ H coefficient function C2,Hg from the usual coefficient function 
C2,g describing the g ^ q transition of light quarks. 
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boundary conditions for evolution. The parton distributions in different quark number regimes 
are related to each other perturbatively, i.e. 

fr'iQ") = E^^'^^^vo ® mQ% (19) 

k 

where the perturbative matrix elements Ajk{Q'^ / m^) jssi containing ln(Q^/m|^) terms are known 
to NNLO, i.e. They relate fiiQ"^) and fi^^iQ'^)-, guaranteeing the correct evolution 

for both regimes. At NLO in the MS scheme they constrain the heavy quarks to begin their 
evolution from a zero value at = m'jj, and the other partons to be continuous at this choice 
of transition point, hence making it the natural choice. 



4.3 General-mass variable flavour number schemes 



The ZM-VFNS has many advantages. However, it has the failing that it simply ignores 0{mfj/Q'^) 
corrections to the coefficient functions, and hence it is inaccurate in the region where is not 
so much greater than mjj. The nomenclature scheme may be thought of as misleading in a 
similar way that zero-mass might. Scheme usually implies an alternative choice in ordering the 
expansion, or a particular separation of contributions between coefficient functions and parton 
distributions, i.e. the inherent ambiguity in a perturbative QCD calculation allows a choice, 
the effects of which become increasingly smaller as higher orders are included. The ZM-VFNS 
misses out 0{m'\j/Q'^) contributions completely, and there is a permanent error of this order. 
This clearly already happens at LO and NLO. The error induced by fitting to HERA structure 
function data using a ZM-VFNS was shown to be up to ~ 6% in the small-x light-quark distri- 
butions by CTEQ in Ref. [8^, resulting in a systematic error of ~ 8% in predictions for vector 
boson production at the LHC This NLO result seems to provide ample evidence for the use of 
a general-mass variable flavour number scheme (GM-VFNS) to provide default parton distribu- 



85 



tions, and indeed this is the approach we have adopted in previous global analyses [12Hl4l. Il8l . |21 
A deflnition of a GM-VFNS was flrst proposed in Re_L 
group have been using the alternative TR scheme 



the ACOT scheme, 



Selas the default since MRST98 [12 



and the MRST 
The 



CTEQ group have since also adopted this convention [SJ, l87 



However, the inherent problems with the ZM-VFNS are thrown into particularly sharp relief 



at NNLO. At 0{aj.) the Ajk 
in the parton distributions. 



are no longer zero at Q 



m 



and lead to discontinuities 



This is similar to the discontinuity at fx 



m 



H 



in asifi"^) at 



NNLO. For the coupling constant this discontinuity is rather small. However, the corresponding 
discontinuities in parton distributions can be signiflcant, about 10% for the gluon, but it turns 
out that c{x,ml) is quite considerably negative at small x, see e.g. Fig. 3 of Ref. [2l|. As it 
happens, the effect of the NNLO massless coefficient function is to make this worse rather than 
better at x ~ 0.0001 — 0.01, and the charm contribution to F2{x, Q"^) is negative at the transition 
point to such a degree that there is an 0(10%) discontinuity in the total structure function at 



as illustrated in Fig. 1 of Ref. 



^When testing our NNLO evolution code against PEGASUS [68| we traced the source of a discrepancy in our 
implementation of Ahq to a typo in Eq.(B.3) of the preprint version of Ref. fssj (fixed in the journal version). 
Specifically, the term Ca Tf (162;+ 16z^) [Li2(— z) + lnzln(l + z)] has f6 rather than 16z in the preprint version. 
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Hence, for a precise analysis of structure functions, and other data, one must use a GM-VFNS 
which smoothly connects the two well-defined limits of < and ^ m'jj. This is already 
the case at lower orders, but becomes imperative at NNLO. However, despite the above reasoning, 
a GM-VFNS is not always used as default even today. This is undoubtedly often because of 
the extra complication compared to the simple ZM-VFNS. However, part of the reason may be 
because the definition is not unique. The reasons for this and the consequences will now be 
discussed, along with a detailed outline of the prescription used in our analysis. A less detailed, 
but more introductory and comparative discussion of the schemes used by MRST/MSTW and 
CTEQ can be found in Ref. 



A GM-VFNS can be defined by demanding equivalence of the nf = n (FFNS) and Uf = n + 1- 
flavour (GM-VFNS) descriptions above the transition point for the new parton distributions 
(they are by definition identical below this point), at all orders, i.e. 

F,(x,Q^) = 5^cjr'"(gvo®A"(Q^) = E^^'"^'wvo®/r'(Q^) (20) 

k j 

The description where the number of active partons is taken to be = n must be identical to 
that when it increases, i.e. nj = n + 1. Hence, the GM-VFNS coefficient functions satisfjQ 

<r'"(QVo = ® ^..(Qvo, (21) 

which, for example, at 0{as) gives for F2{x,Q^): 

c!'r''\QVml) = C^J^T'''\QVml) ® P(J) ln(QVO + <'"^^'^'^(QV0. (22) 

The GM-VFNS coefficient functions, CjJ'^^ , are constrained to tend to the massless limits for 
^ mjj and the Ajk{Q'^ / rri^j) are such that this happens self-consistently. However, the 
cJJ'^^'^^\q'^ /mjj) are only uniquely defined in this massless limit Q'^/m'jj — > oo. For finite 

Q'^/mjj one can swap 0{mjj/Q'^) terms between C^h]5'^^'^°^(<5^/"^h) ^^'^ C'^J'"'^^'^^^(<5^/"^h) 
while maintaining the exact definition in fl22l) . It is clear that this general feature applies to 
all relationships in fl2Tl) . Although the equivalence fl2T]) was first pointed out in general in 
Ref. [s^, and f l22l) is effectively used in defining the original ACOT scheme, the freedom to swap 



0{m'jj/Q'^) terms without violating the definition of a GM-VFNS was first noticed in Ref. 86 
and put to use to define the TR scheme, as described below. This freedom to redistribute 
0{m\/Q'^) terms can be classified as a change in scheme since it leads to an ambiguity in the 
result at a fixed order, but the ambiguity becomes higher order if the order of the calculation 
increases, much like the renormalisation and factorisation scheme (and scale) ambiguities. More- 
over, it is a change of scheme which does not change the definition of the parton distributions. 



^It is implicit that the couphng constant is a function of n flavours on the left-hand side and of n + 1 flavours 
on the right-hand side. 
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only the coefficient functions. This is perhaps a surprising result, which occurs because there is 
a redundancy in fl2T|) . there being one more coefficient function above the transition point than 
below, i.e. that for the heavy quark. 



The original ACOT prescription [85j calculated the coefficient functions for single heavy quark 
scattering from a virtual photon exactly. This might seem the most natural definition. However, 
it assumes that immediately above the transition point a single heavy quark or antiquark can 
exist in isolation. Hence, each coefficient function violates the real physical threshold W"^ > 
4:mjj, since only one heavy quark is produced in the final state, rather than a quark-antiquark 
pair. Moreover, this definition requires the calculation of mass-dependent coefficient functions, 
which becomes progressively more difficult at higher orders. As mentioned above, in the TR 
scheme 86| the ambiguity in the definition of C*J^']l^^'*'°^((5^/m^) was recognised and exploited 



for the ffist time. To be precise, the constraint that {dF^ /dlnQ"^) was continuous at the 
transition point (in the gluon sector) was applied to define the heavy-quark coefficient functions. 
This imposed the correct threshold dependence on all coefficient functions and improved the 
smoothness at = mfj, and did not involve the explicit calculation of new mass- dependent 
diagrams. However, it did involve complications, since it required the convolution of the formal 
inverse of splitting functions with coefficient functions, which itself becomes technically difficult 
at higher orders. 

Since these early definitions there have been various modifications, including a precise def- 
inition of an ACOT-like scheme up to NNLO by Chuvakin, Smith and van Neerven ME A 
major simplification was achieved when the flexibility in the choice of heavy-quark coefficient 



functions was used to define the ACOT(x) prescription [90|, l91j, which in the language used in 
this paper (and in Ref. [s^) would be defined by 



QVO = elz{l + Amyq^) 5(z- ^^f^ ) • (23) 



This gives the LO definition 

(°)(x, Q2) = el (x/x^ax)(^ + ^)(x/x^ax, Q') , (24) 

where Xmax = Q'^/iQ'^ + 4m|^). It automatically reduces to the massless limit C™'j^\z) = 
e^j z 5(1 — z) for Q"^ /m\ oo, and also imposes the true physical threshold 

iy2 _ q2^^ _ y 4^2^_ (25) 



^Indeed in Ref. j89| it is pointed out that at NNLO a further comphcation appears for the first time. There are 
In'^ {Q'^/mD divergences at 0(a|) coming from gluon sphtting into heavy quark-antiquark pairs. These diver- 
gences arise from heavy quark emission diagrams which cancel with opposite ones originating from virtual heavy 
quark loops in the "light quark" coefficient functions. This cancellation is achieved in a physically meaningful 
manner by imposing a cut on the softness of the heavy quark final state in the former process, i.e. the unob- 
servable soft process cancels with the virtual corrections in the "light quark" cross section, and a well-behaved 
remainder enters the "heavy quark" cross section. In practice, after cancellation both contributions are very 
small (both are quark- rather than gluon-initiatcd) at NNLO and we currently include the total in the "light 
quark" sector. In the 0(0^) FFNS they are usually combined instead in the "heavy quark" contribution. 
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This choice of the LO heavy-flavour coefficient function has been adopted in our current pre- 
scription, which we denote the TR' scheme, described in detail in Ref. 72] • For the GM-VFNS 
to remain simple (and physical) at all orders, m, it is necessary to choose 



^2,Hh' {^■•Q'^ l^y) — ^2,Hh' (^/^^max), (26) 

which is the implicit ACOT(x) definition, and is our choice. It removes one of the sources of 
ambiguity in defining a GM-VFNS. However, there are others. 

One major issue in a complete definition of the GM-VFNS, is that of the ordering of the 
perturbative expansion. This ambiguity comes about because the ordering in as for F^(x,(5^) 
is different for the number of active flavours nj = n and Uf = n + 1 regions: 

n — flavour n + 1 — flavour 



NNLo )^ c!^:^^ ® /; {^ye, c^sr''^'^ ® fr' 

with obvious generalisation to even higher orders. This means that switching directly from an n 
flavours fixed order to the n + 1 fixed order leads to a discontinuity in F^{x, Q"^). As with the 
discontinuities in the ZM-VFNS already discussed this is not just a problem in principle — the 
discontinuity is comparable to the errors on data, particularly at small x. 

Hence, any definition of a GM-VFNS must make some decision how to deal with this, and 
the ACOT-type schemes have always made a different choice to that for the TR-type schemes 
used in our analyses. The ACOT-type schemes simply define the same order of 0:5 both below 
and above the transition point. For example at NLO the definition is 

F«(.,<3^) = ^Cf-'"«9" - g (C,™* ® (ff + 77) +C™' »<,"«). (27) 

This clearly maintains continuity in the structure function across the transition point. However, 
it only contains information on LO heavy flavour evolution below = m|^, since C^^^^'^^^ 
only contains information on the LO splitting function, but the heavy quarks evolve using NLO 
split ting functions above = — a big change at small x. The TR scheme, defined in 
Ref. 86|, and all subsequent variations used in our analyses, try to maintain the correct ordering 



in each region as closely as possible. For example at LO we have the definition 

Ff(x,Q^) = ^^C!X'^'\QVml)(^g-{Q^) (28) 

4:71 ' ^ 

i.e. we freeze the 0{as) term when going upwards through = m?^. This generalises to higher 
orders by freezing the term with the highest power of q;^ in the definition for < m\ when 



18 



moving upwards above m|^. Hence, the definition of tlie ordering is consistent witliin eacli region, 
except for the addition of a constant term (which does not affect evolution) above = mj^, 
which becomes progressively less important at higher Q^, and whose power of as increases as 
the order of the perturbative expansion increases. 



we 



This definition of the ordering means that in order to define our VFNS at NNLO [72 
need to use the heavy-flavour coefficient functions for < m\ (and this contribution 

will be frozen for > rnjj). As mentioned above, these coefficient functions are not yet 



calculated. However, as explained in Ref. [72], we can model this contribution using the known 
leading threshold logarithms [o^] and leading ln(l/a;) terms derived from the /cr-dependent 
impact factors [i^. This results in a significant contribution at small and x with some 
model dependence. However, variation in the free parameters does not lead to a large change, 
as discussed in Section [9]0 Up to this small model dependence we have a full NNLO GM-VFNS 
with automatic continuity of structure functions across heavy flavour transition pointsl^ This 
is certainly the most complete treatment of heavy-flavour effects currently used in any NNLO 
analysis. 

Scheme dependence 

Although all definitions of the GM-VFNS become very similar at very high Q^, the difference in 
choice can be phenomenologically important. For example, our definition effectively includes ex- 
actly one higher order than ACOT-type schemes for < rnjj, and the value of this contribution 
at = m\ is carried to higher Q^. Since at small x and near threshold the higher orders in 05 
are accompanied by large corrections, this leads to large differences below the transition point, 
which are still irapqrtant a little way above the transition point. This is shown for -F|'^(x, Q'^) 



in Fig. 2 of Ref. [72| where the two choices are shown for the same parton distributions. More 
clearly, this difference in the definition of the ordering is the main difference in the NLO predic- 
tions from MRST and CTEQ in the comparison to HI data on F^{x^ Q'^) jisj, shown in Fig. 4 
of the same paper [72]. 



The inclusion of the complete GM-VFNS in a global fit at NNLO first appeared in Ref. [21 
and led to some important changes compared to our previous NNLO analyses, which had a much 
more approximate inclusion of heavy flavours (which was explained clearly in the Appendix 



of (l5|). A consequence of including the positive C(a|) coefficient functions at low is that 
the NNLO F|'^(x, Q^) automatically starts from a higher value at low Q^. However, at high 
Q^, the structure function is dominated by (c + c)(x, Q^)- This has started evolving from a 



^It should be stressed that this model is only valid for the region < rnj^, and would not be useful for 
a full NNLO FFNS since it contains no information on the large Q^/mj^ limits of the coefficient functions. 
A more general approximation to the 0{ag) coefficient functions could be attempted, but full details would 
require first the calculation of the 0(a|) matrix element Ahq- This more tractable project is being investigated 
at present 94 1. An approximation to the a| coefficient functions using logarithmically-enhanced terms near 



threshold and exact scale-dependent terms has recently been proposed in Ref. ^95|] . 

^There are actually 0(a|) discontinuities due to terms such as C^^^^^'^^^ i^{H + H) and C^^/'g ^^'^^^ ® 5""*""^, 
i.e. 0{as) coefficient functions convoluted with 0(q;|) discontinuities in partons. These would be cancelled at 
NNNLO by discontinuities in 0{a^) coefficient functions. In practice the imposition of the correct threshold 
behaviour in all coefficient functions minimises these effects and they are very small. 
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significantly negative value at = m^. The parton distributions in an NNLO fit readjust so 
that the light flavours evolve similarly to those at NLO, in order to fit the data. Since the 
heavy flavour quarks evolve at the same rate as light quarks, but at NNLO start from a negative 
starting value, they remain lower than at NLO for higher Q^. Hence, there is a general trend: 



Ff{x,Q'^) is flatter in at NNLO than at NLO, as shown in Fig. 4 of Ref. [2l|. It is also 
flatter than in our earlier (approximate) NNLO analyses. We found that this had an important 
effect on the gluon distribution when we updated our NNLO analysis. As seen in Fig. 5 of 



Ref. [2l|, it led to a larger gluon for x ~ 0.0001-0.01, as well as a larger value of a5(M|), both 
compensating for the naturally flatter evolution, and consequently leading to more evolution of 
the light quark sea. Both the gluon and the light quark sea were up to 6-7% greater than in the 



2004 set 181] for = 10"^ GeV^, the increase maximising at a; = 0.0001-0.001. As a result there 
was a 6% increase in the predictions for aw and az at the LHC. This surprisingly large change 
is a correction rather than a reflection of the uncertainty due to the freedom in choosing heavy 
flavour schemes. The treatment of heavy flavour at NNLO is the same in the sets presented 
in this paper as for the MRST 2006 set, and as we will see, there are no further changes in 
predictions of the same size as this increase in going from the 2004 to the 2006 analyses. This 
demonstrates that the MRST 2004 NNLO distributions should now be considered to be obsolete. 



Our 2006 NNLO parton update [21j was made because this was the first time the heavy 
flavour prescription had been treated precisely at NNLO and also because there was previously 
no MRST NNLO set with uncertainties. The data used in the analysis were very similar to the 
2004 set, and since a consistent GM-VFNS was already used at NLO, and a set with uncertainties 
already existed, no new corresponding release of a NLO set was made along with the 2006 NNLO 
set. With the benefit of hindsight, it is interesting to check the effect on the distributions due 
to the change in the prescription for the GM-VFNS at NLO without complicating the issue by 
also changing many other things in the analysis. To this end we have obtained an unpublished 
MRST 2006 NLO set, which is fit to exactly the same data as the MRST 2006 NNLO set@ 

The comparison of the up quark and gluon distributions for the MRST 2006 NLO set and 



the MRST 2004 NLO set, i.e. the comparable plot to Fig. 5 of Ref. [2l| for NNLO, is shown in 
Fig. [3l As can be seen it leads to the same trend for the parton distributions as at NNLO, i.e. an 
increase in the small-x gluon and light quarks, but the effect is much smaller — a maximum of 
a 2% change. Also, the value of the coupling constant increases by 0.001 from the 2004 value of 
asiM"^) = 0.120. Again, this is similar to, but smaller than, the change at NNLO. Hence, we 
can conclude that the change in our choice of the heavy-flavour coefficient function alone leads 
to changes in the distributions of up to 2%, and since the change is simply a freedom we have in 
making a definition, this is a theoretical uncertainty on the parton distributions, much like the 
frequently invoked scale uncertainty. Like the latter, it should decrease as we go to higher orders. 
The ambiguity simultaneously moves to higher order, but it is difficult to check this explicitly 
since our main reason for making our change in the choice of heavy-quark coefficient functions 



was the difficulty of applying the original procedure in Ref. [86|] at NNLO. Certainly an absolute 



maximum of 2% of the 6-7% change, in the predictions for aw and az at the LHC in going from 



^We do not intend to officially release the MRST 2006 NLO set, since it is superseded by the present MSTW 
2008 NLO analysis, but it is nevertheless available on request. 
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(a) 



Gluon distribution at = 10" GeV^ 



MRST 2004 NLO (90% C.L.) 
MRST 2006 NLO (unpublished) / 





Figure 3: A comparison of the unpublished MRST 2006 NLO parton distributions to the MRST 
2004 NLO distributions at = 10^ GeV^. In order to illustrate the significance of the size of 
the differences, the fractional uncertainty on the MRST 2001 distributions is used for the 2004 
distributions. The corresponding comparison at NNLO can be seen in Fig. 5 of Ref. 21 . 



the 2004 to the 2006 NNLO parton sets, is due to true ambiguities, and the remaining 5% is due 
to the correction of the flaws in the previous approach. 

To close this section, we note that the MRST 2006 NLO distributions lead to a nearly 3% 
increase in the predictions for aw and az at the LHC compared to MRST 2004 NLO, though 
there is very little change at the Tevatron, where the typical values of x probed are nearly an 
order of magnitude higher. As with the distributions themselves, this variation is a similar size to 
the quoted uncertainty in the cross sections, and again this is a genuine theoretical uncertainty. 
Our most up-to-date predictions for LHC and Tevatron W and Z total cross sections, and 
comparison to the various previous sets, will be made in Section [T5l 



Longitudinal structure function 

One also has to be careful in defining the GM-VFNS for the longitudinal structure function. In 
one sense this is not so important since the contribution from the longitudinal structure function 
to the total reduced cross section measured in DIS experiments is small, and the errors on direct 
measurements of Fl are comparatively large. However, the importance is increased by the larger 
ambiguity inherent in the definition of Fl as compared to F2. 

This large ambiguity occurs because if one calculates the coefficient function for a single 
massive quark scattering off a virtual photon, there is an explicit zeroth-order contribution 

-el. ^) . (29) 

This disappears at high and the correct zero-mass limit is reached. 

Such a zeroth-order coefficient function is implicit in the original ACOT definition of a GM- 
VFNS, and leads to a peculiar behaviour of Fl'^ just above = m^. It is convoluted with the 
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heavy flavour distribution, which for Q"^ just above is small in magnitude. However, the 
coefficient function is large near m^, while the unsubtracted (i.e. FFNS) gluon and singlet-quark 
coefficient functions are suppressed by a factor of v^, where v is the velocity of the heavy quark 
in the centre-of-mass frame, and are very small for low Q^. This means that this zeroth-order 
heavy-flavour contribution dominates just above = m^, despite the fact that in a FFNS, 
where the cc pair has to be created, as it must in reality, the contribution is absent. 

The contribution from cj^^"^ ® (c + c) turns on rapidly just above m^, dominating other 
contributions, then dies away as mj^/Q"^ becomes small. This leads to a distinct bump in 
F^^{x,Q'^) for Q'^ just above ml, as pointed out in Ref. 8^. In principle this cancels between 



orders in a properly defined GM-VFNS, as this contribution implicitly appears in the subtraction 
terms for the gluon and singlet-quark coefficient functions with opposite sign to its explicit 
contribution. However, the cancellation is imperfect at finite order, and even the partially 
cancelled contribution dominates at NLO. If this coefficient function is implemented it leads to 
peculiar behaviour for slightly above m^. At NNLO, where heavy-flavour distributions begin 
at = ml with negative values, the "bump" in F^{x, Q^) is negative, as illustrated in Fig. 18 
of Ref. [89|, highlighting the unphysical nature of this contribution. 

Hence, as in Ref. [86j we choose to ignore the explicit single heavy-quark-photon scattering 
results. We define the longitudinal sector in what seems to us to be the most physical gener- 



alisation of the definition for F^lxyQ"^), as explained in Ref. 72]. The heavy-quark coefficient 
functions are simply those for the light quarks, with the upper limit of integration moved from 
1 to Xmax = Q^/{Q^ + 4m|^). Thus the physical threshold of > 4m|^ is contained in all 
terms, and there are no spurious zeroth-order terms. These could only make a contribution 
if one works in the framework of single heavy quark scattering in the region of low where 
the parton model for the heavy quark is least appropriate. The definition of the SACOT-type 



scheme [96|, and particularly the SACOT(x) scheme used in the global fits [84|, also avoids this 
undesirable zeroth-order coefficient function. We will discuss our predictions for Fl{x,Q'^) in 
Section I 



Charged-current structure functions 

The extension of the GM-VFNS to charged currents is most important for the heavy-flavour 
contribution to neutrino structure function data from CCFR, NuTeV and CHORUS, where the 
dominant process is z/ + s — > £~ -|- c (with small Cabibbo mixing) and the charge conjugate 
process. As such, this is particularly important for an analysis of the dimuon production data 
from CCFR and NuTeV, discussed in more detail in Section [HI There are also charged-current 
data from HERA, but this is less precise (though it does not need nuclear corrections, unlike 
the neutrino data), and is at sufficiently high that the heavy quarks are effectively massless, 
though the GM-VFNS is applied in this case also. 

The general procedure for the GM-VFNS for charged-current deep-inelastic scattering works 
on the same principles as for neutral currents — one can now produce a single charm quark 
from a strange quark so the threshold is now at Xmax = Q"^ / iQ"^ + ^c)- However, as explained 
in Ref. [72|, there is a complication because the massive FFNS coefficient functions are not 



known at 0{a^g) (only asymptotic limits (97| have been calculated). These coefficient functions 
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Figure 4: The gluon coefScient function at 0{as) in NC and CC processes. 



are needed in our GM-VFNS at low Q"^ at NLO, and at all at NNLO — though in the 
latter case the definition of the GM-VFNS means that the ln((5^/m^) terms are subtracted, and 
the 0{m1/Q'^) terms die away at high Q^, so the GM-VFNS coefficient functions tend to the 
precisely known massless limits for large Q^/m^. 

The initial proposal to deal with this, outlined in Ref. [t^], was to assume that the mass- 
dependence in the 0{a\) coefficient functions is the same as for the neutral-current functions, 
but with the threshold in Ami replaced by a threshold in m^. It was noted that this meant that 
the coefficient functions at least satisfy the threshold requirements, and tend smoothly to the 
correct massless limits, so were very likely to be an improvement on the ZM-VFNS. However, in 
the course of the analyses performed in the latest global fit we have noticed various complications. 
One consideration is that the neutrino cross sections are given by expressions of the form 

^ « (i-y + ^-^)i^2""(a:,Q^)-fFf^(x,g^)±y(i-|)xF3c^(x,g^), (30) 

where y ~ 0.3 — 0.8. Unlike the general case for neutral currents, where the F2{x, Q"^) term dom- 
inates, all three terms are important for charged currents and in the heavy-flavour contribution 
there are signiflcant cancellations between them, so all need to be treated very carefully. 

Let us flrst discuss F2^'"(x, Q^). At O (a'i ) this is dominated by the gluon contribution. 



The simple prescription, suggested in Ref. [72|, gives large corrections to F2^{x,Q'^) at low Q"^ 
because the lower threshold compared to the neutral-current case leads to a longer convolution 
length. However, let us consider the comparison of the gluon coefficient functions at 0{as), 
represented by the two diagrams in Fig. HI 

The neutral-current coefficient function is infrared finite and positive. On the other hand the 
charged-current coefficient function diverges due to the coUinear emission of a light (i.e. strange) 
quark. After subtraction of this divergence, via the usual factorisation theorem in the MS 
scheme, there is an approximate factor of (1 + 21n(l — z)) in the coefficient function at low Q^. 
It is negative for high z even in the FFNS. At higher Q^, the finite parts of the neutral-current 
and charged-current GM-VFNS coefficient functions, after the subtraction of the ln((5^/m^) 
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terms, converge to the same massless expression, but from qualitatively different forms at lower 

At the neutral-current FFNS coefficient function is positive with threshold log en- 

hancements. On the other hand, obtaining the approximate charged- current expression by just 
a change in threshold limits leads to a large contribution, which, bearing in mind the above 
discussion, is unlikely to be accurate. At (9(a|), the charged- current diagram will still have an 
extra emission of a strange quark and a corresponding coUinear subtraction. So we choose to 
obtain the charged-current coefficient function by both the change in threshold kinematics and 
by introducing a (1 + 2(m^/(5^) ln(l — z)) factor. This still leads to the correct large Q"^ /ml 
limit, but is very likely to be a more accurate representation of the low region. The same 
procedure is applied to the similar singlet-quark contribution. 

We now consider xF^'^{x,Q'^) at 0{a'g). For massless quarks this contribution is zero for 
initial gluons and singlet quarks, but the coefficient function C^^'"-'^(a;, m^, Q^) is non-zero for 
finite Q'^ /ml. However, it must vanish as both Q"^ /ml oo and W'^/ml oo. Hence, our 
model for C^^'^^'^'^\x,ml,Q'^) is weighted by a factor ml/W"^ = mlz/({l — z)Q'^), as is the 
singlet-quark contribution. It is important that a suppression of this type is implemented. 
Otherwise the contribution is potentially anomalously important at low and x. 

Finally, there is a complication in the ordering for the longitudinal charged-current heavy 
fiavour production. In the massless limit the lowest-order contribution, simplified by neglecting 
the Cabibbo-suppressed rf-quark contribution, is 

F^^'^^x, Q') = 2 ® g{x, Q') + ^^(a;) ® s{x, Q')). (31) 

However, for a massive quark there is a zeroth-order contribution 

Ff^'^(x,Q^) = 2-^e^(e,Q^), (32) 

where ^ = x(l -|- ml/Q"^). Note that this is unlike the neutral-current case for Fl{x,Q^), 
where there was also a zeroth-order contribution. Here it is due to a real physical process, 
i.e. + s c, rather than one which only makes sense in the limit where the charm quark 
is most definitely behaving like a massless parton. Hence, for the charged-current case, the 
zero-order contribution must be included. This means that there is a difference in orders below 
and above the transition point, i.e. the FFNS begins at zeroth order whereas the ZM-VFNS 
begins at first order — opposite to the case for the neutral-current F2{x,Q^). Again a choice 
in ordering must be made. We choose to obtain the correct limits in both regimes along with 
maintaining continuity. That is, we use (1521) to define the LO contribution for < ml, whereas 
for > ml the LO contribution is defined by 

^'^(x, Q') = 2 e^(e, Q') + 2 (^1 - asiCU^) ® gix, Q') + C^,,(x) ® six, Q')). 

(33) 

At high the first term naturally dies away leading to the normal massless limit. We can easily 
generalise the prescription to higher orders by including the next term in the as expansion on 



24 



both sides of the transition point, with the (1 — ml/Q^) factor always multiplying the highest- 
order term in the region above the transition point. 

In principle, we could also make some use of the 0{ag) charged-current coefficient functions 
in the same way as we do for neutral currents. However, the region where they make a significant 
impact is overwhelmingly in the small x and regime accessed only by HERA neutral-current 
measurements. The amount of modelling required for these terms in charged-current processes 
is large. Since they are unlikely to have much effect at all, we simply omit them. 

With the chosen prescription for charged-current heavy fiavour described above, we find that 
the description of the data, and the resulting parton distributions, are rather stable in going 
from LO — > NLO NNLO. This will be discussed further in later sections. 



4.4 Intrinsic heavy flavour 

Throughout the above discussion we have made the assumption that all heavy quark (c, b) fiavour 
is perturbatively generated from the gluon and lighter fiavours. Indeed, this has always been 
the assumption in our analyses. It is justified up to corrections of 0{AQ^j^/m'jj). This potential 
power-suppressed correction is known as intrinsic heavy fiavour. Although it is parametrically 
small it may be relatively enhanced at high values of x 98[, and a plausible measure of the 



contribution is that at low scales the integrated number density to find intrinsic charm is < 1%. 
There is little within the current global fit which can tightly constrain this possible intrinsic 
fiavour contribution. The possibility of different types of intrinsic charm contribution has been 
studied recently by CTEQ [99]. They conclude that a global fit can allow up to a factor ~ 3 
times more high-x enhanced intrinsic charm than the traditional value of 1%, and that a sea-like 
contribution can carry a momentum fraction up to 2% of that of the proton (the integrated 
number density being infinite in this case). Indeed, for the sea-like distribution, they find a 
preference for a non-zero value, but we believe that this merely compensates for the terms which 
are systematically absent in the ACOT(x) scheme at low compared to our TR' scheme or 
to the FFNS. There seems to be no theoretical impetus to consider a sea-like intrinsic heavy 
flavour to be anything other than tiny. We will consider the possibilities for a high-x enhanced 
contribution in Section [91 However, in order to do so, we discuss briefly how it flts into the 
GM-VFNS framework. 

Intrinsic flavour and GM-VFNS deflnitions were discussed in Ref. [Tool . We briefly recall 



them here. Allowing an intrinsic heavy quark distribution actually removes the redundancy in 
the deflnition of the coefficient functions in the GM-VFNS, and two different deflnitions of a 
GM-VFNS will no longer be identical if formally summed to all orders, though they will only 
differ by contributions depending on the intrinsic flavour. The reason is as follows. Consider 
using identical parton distributions, including the intrinsic heavy quarks, in two different flavour 
schemes. The heavy-quark coefficient functions at each order are different by 0{m\j/ Q"^). This 
difference has been constructed to disappear at all orders when combining the parton distribu- 
tions other than the intrinsic heavy quarks, but will persist for the intrinsic contribution. The 
intrinsic heavy-flavour distributions are of (9(AQ^^/m|^), and when combined with the differ- 
ence in coefficient functions the mass-dependence cancels leading to a difference in structure 



functions of 0{h.Q(jj-)/Q ). It has been shown |101| that for a given GM-VFNS the calculation of 
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the structure functions is limited in accuracy to 0{AQf^<j^/Q^). Hence, when including intrinsic 
charm, the scheme ambiguity is of the same order as the best possible accuracy one can obtain 
in leading-twist QCD, which is admittedly better than that obtained from ignoring the intrinsic 
heavy flavour (if it exists) as increases above mjj. It is intuitively obvious that better accu- 
racy will be obtained from defining a GM-VFNS where all coefficient functions respect threshold 
kinematics, e.g. the ACOT(x) and TR' schemes, than from some older schemes which do not, 
and which violate physical thresholds when combined with intrinsic heavy-flavour contributions. 
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5 Global parton analyses 



In this section we present a summary of three global PDF analyses performed using the theo- 
retical formalism outlined in Section [31 together with the treatment of heavy flavours given in 
Section HI We perform fits at LO, NLO and NNLO starting from input parton distributions 
parameterised at Qq = 1 GeV^ in the form shown in Eqs. fl6|)- f|T2|) . 



5.1 Choice of data sets 



The data sets included in the analyses are listed in Table [21 together with the values, defined 
corresponding to each individual data set for each of the three fits, and the 



in Section 15.2 

number of individual data points fitted for each data set. The data sets in Table [2] are ordered 
according to the type of process. First we have the fixed-target data, which are subdivided into 
I'^N structure functions, low-mass Drell-Yan (DY) cross sections, and vN structure functions 
and dimuon cross sections. Then we list data collected at HERA, and finally data collected at 
the Tevatron. More detailed discussion of the description of the individual data sets will be 
given later in Sections [714T21 

For the DIS data we make cuts on the exchanged boson virtuality oi > 2 GeV^ and on 
the squared boson-nucleon centre-of-mass energy, W"^ = Q^{l/x — 1) +m%, of W"^ > 15 GeV^ at 
NLO and NNLO, while W"^ > 20 GeV^ at LO due to the absence of higher-order large-x terms. 
The cut on W"^ is necessary to exclude data for which potentially large higher-twist contributions 
are present, such as target mass corrections [115j . For the v{i')N xF^ data we only include data 



with W"^ > 25 GeV^ for all three fits. The higher-twist contributions for F2 {x, Q^) do not 
contribute to the Adler sum rule, i.e. /o^dx F^'^(x,Q'^) = 0, which means they must be well 
behaved as x 0. However, there is no such restriction on the hi gher-twist corrections for xF-- 



and renormalon calculations imply that they are large jll6 . 117]. The CHORUS data extend 



into the region of expected large higher-twist corrections, so the above c ut is nece ssary to avoid 



contamination (which is qualitatively similar to that predicted in Refs. |116 
As i n pr evious MRST fits, data on the deuteron str uctu re function 
F2/F2 |l03 |. are corrected for nuclear shadowi ng e ffects I1I8 



na). 



and the ratio 
No such corrections are ap- 
plied to the E866/NuSea pd/pp Drell-Yan ratio [lQ8j, though we expect such corrections to be 



small compared to the experimental errors, and it is not clear how such corrections would be 
applied. For other nuclear target data we apply an improved nuclear correction discussed in 
detail in Section I7.3[ We apply a correction of —3.4 % to the luminosity of the published HI MB 



97 data [l09l | following a luminosity reanalysis [119 



The Tevatron high-pr jet data are included in the fit at each order. The full NNLO co rrect ions 



are not known in this case, but the approximation based on threshold corrections 120l 



included in the fastnlo package 75 



are 

There is no guarantee that these give a very good 
approximation to the full NNLO corrections, but in this case the NLO corrections themselves 
are of the same order as the systematic uncertainties on the data. The threshold corrections are 
the only realistic source of large NNLO corrections, so the fact that they provide a correction 
which is smooth in px and small compared to NLO (and compared to the systematic uncertainties 
on the data) strongly implies that the full NNLO corrections would lead to very little change 
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Data set 


LO 


NLO 


NNLO 


BCDMS HP F2 [32] 


165 / 153 


182 / 163 


170 / 163 


BCDMS ^id F2 [102] 


162 / 142 


190 / 151 


188 / 151 


NMC HP F2 [33] 


137 / 115 


121 / 123 


115 / 123 


NMC lid F2 [33] 


120 / 115 


102 / 123 


93 / 123 


NMC iJ.n/iip [1031 


131 / 137 


130 / 148 


135 / 148 


E665 HP F2 [104] 


59 / 53 


57 / 53 


63 / 53 


E665 /id F2 [104] 


49 / 53 


53 / 53 


63 / 53 


SLAG ep F2 [105, 106] 


24 / 18 


30 / 37 


31 / 37 


SLAG ed F2 [105, 1061 


12 / 18 


30 / 38 


26 / 38 


NMG/BGDMS/SLAG Fl [32-34] 


28 / 24 


38 / 31 


32 / 31 


E866/NuSea DY [107] 


239 / 184 


228 / 184 


237 / 184 


E866/NuSea pd/pp DY [108] 


14 / 15 


14 / 15 


14 / 15 


NuTeV uN F2 [37] 


49 / 49 


49 / 53 


46 / 53 


GHORUS vN F2 \38\ 


21 / 37 


26 / 42 


29 / 42 


NuTeV i/TV XF3 [32 


62 / 45 


40 / 45 


34 / 45 


GHORUS i/A^ xFg [38] 


44 / 33 


31 / 33 


26 / 33 


GGFR j/A -> /i/^A [39] 


63 / 86 


66 / 86 


69 / 86 


NuTeV i/A ^ /i/iA [39] 


44 / 40 


39 / 40 


45 / 40 


HI MB 99 e+p NG [31] 


9/8 


9/8 


7/8 


HI AIB 97 e+p NG [109] 


46 / 64 


42 / 64 


51 / 64 


HI low Q2 96-97 e+p NG [^09] 


54 / 80 


44 / 80 


45 / 80 


HI high g2 98-99 e-p NG [110] 


134 / 126 


122 / 126 


124 / 126 


HI high g2 99-00 e+p NG [35] 


153 / 147 


131 / 147 


133 / 147 


ZEUS SVX 95 e+p NG [111] 


35 / 30 


35 / 30 


35 / 30 


ZEUS 96-97 e+p NG [112] 


118 / 144 


86 / 144 


86 / 144 


ZEUS 98-99 e-p NG [113] 


61 / 92 


54 / 92 


54 / 92 


ZEUS 99-00 e+p NG [114] 


75 / 90 


63 / 90 


65 / 90 


HI 99-00 e+p GG [35] 


28 / 28 


29 / 28 


29 / 28 


ZEUS 99-00 e+p GG [36] 


36 / 30 


38 / 30 


37 / 30 


Hl/ZEUS ep i^2=harm [41_47] 


110 / 83 


107 / 83 


95 / 83 


HI 99-00 e+p inch jets [59] 


109 / 24 


19 / 24 


— 


ZEUS 96-97 e+p inch jets [57] 


88 / 30 


30 / 30 




Zj£l/Uo yo— uu e p mci. jeijs 


1 no / "in 


1 7 / "^n 

1 ( / oU 




D0 H pp incl. jets [56| 


193 / 110 


114 / 110 


123 / 110 


GDF H pp incl. jets \5£ 


143 / 76 


56 / 76 


54 / 76 


GDF IIW ^Iv asym. [48] 


50 / 22 


29 / 22 


30 / 22 


D0 H W ^ asym. [49] 


23 / 10 


25 / 10 


25 / 10 


D0 H Z rap. [53] 


25 / 28 


19 / 28 


17 / 28 


GDF H Z rap. [52] 


52 / 29 


49 / 29 


50 / 29 


All data sets 


3066 / 2598 


2543 / 2699 


2480 / 2615 



Table 2: The values of x^/^pts. for the data sets included in the global fits. For the NuTeV uN — > 
fifiX data, the effective number of degrees of freedom is quoted instead of A^pts. since smearing 
effects mean that nearby data points are highly correlated |39j. The details of corrections to 
data, kinematic cuts applied and definitions of s'^re contained in the text. 
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in the PDFs. Since these jet data are the only good direct constraint on the high-x gluon we 
choose to include them in the NNLO fit judging that the impact of leaving them out would be 
more detrimental than any inaccuracies in including them without knowing the full NNLO hard 
cross section. The HERA jet data, for which the NNLO coefficient functions are also unknown 
at present, are omitted from our NNLO fit. This is because they do not provide as important a 
constraint as the Tevatron jet data, and also because the NLO corrections are much larger than 
for Tevatron jet data (due to the lower scales) and there is not even some partial information on 
the NNLO result available. However, we compare the HERA jet data to results calculated with 
the NNLO PDFs using the NLO coefficient functions in Section [121 finding surprisingly good 
agreement. 

5.2 definition 

The global goodness-of-fit measure is defined as Xgiobai ~ J2nXn^ where, in all previous MRST 
fits, 

S/T 1 \r \ — \^ iDn,i — Tn,i{{a})/Afn) 2 /qa\ 

X„lia|,yV„j - 2^ / uncorr.p _^ fa™":)^ 

Here, n labels a particular data set, or a combination of data sets, with a common (fitted) 
normalisation A/'„, i labels the number of individual data points in that data set, and k labels 
the individual correlated systematic errors for a particular data set. The individual data points 
Dn^i have uncorrelated (statistical and systematic) errors crJJ'^-'^™''- and correlated systematic errors 
^riki- The theory prediction T„ j({a}) depends on the input PDF parameters {a}. The MRST 
fits added all uncorrelated and correlated experimental errors in quadrature in the definition 
except for Tevatron jet production. Here, the full correlated error information was used 
but the procedure was more complicated; see Section [121 The second term in (IMj) is a possible 
penalty term to account for deviations of the fitted data set normalisations Mn from the nominal 
values of 1. (In practice, the MRST fits took = 0, although constraints were applied as 
explained below.) We will now discuss improvements to the treatment of data set normalisations 
made in the present analysis, then we will discuss improvements in the treatment of correlated 
systematic errors. 

Treatment of data set normalisations 

In the MRST analyses, the normalisations of several data sets were fioated and then fixed, 
checking that these normalisations lay inside the quoted one-sigma error range. For example, 
the HI data sets from different running periods were all taken to have a common normalisation 
of exactly 1, while the ZEUS data sets from different running periods were all taken to have a 
common normalisation of 0.98. In general, HERA data from different runs are allowed to have 
different normalisations, therefore we now split the HI and ZEUS data into different running 
periods, and allow a different normalisation for each of these running periods^ The normalisa- 

^°We take the normalisations of the HI and ZEUS data on F-f^ to be fixed at f , although here the experimental 
errors are anyway much larger than for the total reduced cross sections, and the Fl'^ data do not provide a strong 
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tions of the data sets are now all allowed to go free at the same time as the PDF parameters, 
including in the determination of the PDF uncertainties discussed in Section [61 

Note that the (fitted) normalisation parameters Mn are included in (IMj) in the correct manner 
to avoid bias in the result, i.e. if the central value of the data is rescaled by a factor Mn then 
so are the uncertainties. If fitting a single data set with no normalisation constraint on the 
theory we would automatically obtain Mn = 1- However, the fit does allow data sets to have 
different relative normalisations, and additionally the sum rules on the PDFs do influence the 
normalisation to some extent. 

The MRST analyses did not include any penalty term in the definition for data set 
normalisations Mn differing from the nominal values of 1. The usual choice for the penalty term 
for the data set normalisations is 

(35) 

where cr^ is the one-sigma normalisation error for data set n. Alternatively, it has been proposed 
that normalisati on u ncertainties should behave according to a log-normal distribution, where the 



penalty term is |121 



Xl = In (Mn) \ 3 + I , (36) 

which reduces to the usual quadratic term (15^ for small and for small deviations of Mn 
from 1. We find that, using fl35|) . the best-fit data set normalisations tend to stray outside 
their nominal one-sigma range, with all the largest shifts being in the downwards direction. 
It has long been known that both LO and NLO fits woul d pr efer to have more than 100% 



momentum for the PDFs if allowed (e.g. at NLO in Ref. 122], the extra momentum in the 



fit regions for the conservative NLO PDFs in Ref. [17|, and the large momentum violation 



at LO seen in Ref. |66l]) so imposing momentum conservation on PDFs conversely leads to a 
preference for lower data set normalisations. However, this is a diminishing effect as we increase 
the order of the QCD calculation, implying that it is an artifact of an incomplete theory at 
lower orders, and as such is an undesirable systematic effect. Additionally, it has been claimed 
that normalisation uncer tainties are expected to behave more like a box shape than the usual 



Gaussian behaviour |123l | , and indeed, the term in does alter the shape of the penalty in the 
downwards direction in this type of manner, albeit to a very small degree for errors of only a few 
percent. Hence, taking into account the systematic effect at low perturbative orders and that 
normalisations are very unlikely to be perfectly Gaussian, we use a more severe quartic penalty 
term for the normalisations, i.e. 

2 f l-Mn Y for, 

This binds the normalisations more strongly to the range Mn £ [1 — cr^, 1 + c"^], but in practice 
it is far more the case that it stops them from floating too low. If more than one data set has 
a common normalisation, the penalty term is divided amongst those data sets according to 
the number of data points. 



constraint on the PDFs. 



30 



The resulting (re) normalisation factors Mn for each fit are listed in Table [31 



Correlated systematic errors 



The insensitivity of both the central values of the PDFs and their uncertainties to the complete 
inclusion of correlation information for the DIS data was confirmed in the benchmark PDF fits 
in the proceedings of the HERA-LHC workshop 



71 



and an extensive discussion of the effect 

To summarise the 



14 



of correlations for HERA data was presented in the appendix of Ref. 
latter, the main effect was by far the systematic shifts due to changes in normalisation which 
we now treat in full (ind eed a strong hint of the necessity for the normalisation correction 119 



of the HI MB 97 data [io9| was seen). Beyond this, the introduction of full correlations led 
to rather small changes and those were, in general, precisely such as to flatten the evolution of 
F2{x, Q^) allowing the gluon to reduce. As previous comments on normalisations and momentum 
conservation suggest, this is itself likely correlated to shortcomings of the theory at low orders. 
Hence, we maintain that the lack of the correlations has little effect on our results, and continue 
to make this simplification. This i s par ticularly the case since we note that for the preliminary 
averaged HERA cross section data |l24 |, sources of correlated uncertainty are often dramatically 
reduced (sometimes by a factor of 3-4), and in the averaging between HI and ZEUS the data 
points frequently move relative to each other in a rather different manner than a data set does 
relative to theory in a fit. We intend to include full correlations for the averaged HERA data, 
when they are published and will be better understood, but at this stage they will also be very 
small indeed. We also note that in some cases a textbook treatment of correlated uncertainties 



can lead to peculiar results. Indeed in Ref. [125[ | it was seen that the high-y turnover in reduced 
cross section data due to the influence of Fl{x,Q'^) could be eliminated completely in a fit by 
letting the correlated systematic uncertainty due to the photoproduction background move by 
about two-sigma. This seems unlikely (though not impossible), and it is certainly the case that 
the distribution of uncertainties in this background are far from Gaussian, so the conventional 
treatment may give surprising results. It is definitely the case that the absolute values for 
the best fit do change depending on how the correlations are treated, but given that our new 
tolerance determination, described in detail in Section [^21 relies on changes in relative to the 
best fit, we are confident that we are even less sensitive to such details than in previous studies. 



e.g. the aforementioned [ly, |71 



Nevertheless, for selected newly-added d ata s ets we do include the full correlated error infor- 
mation. Instead of (|34ll . the is given by 126 



Xl{{a},^fn) = Yl 



=1 



' Dn,^-TnA{a})/^fn \ 
uncorr. I 



+ 



E 



n,k 



+ 



(38) 



where Dn,i = -D„,i - X)fc=i 
error s in order to give the best fit 



corr. 

' n,k '-'n.k.i 



that 126 



are the data points allowed to shift by the systematic 
Minimising Xn with respect to rn,k gives the analytic result 



r-„,fc({a},7V;) = ^(A ^)kk'Bk'{{a},K.] 



(39) 
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1 


1 


1 


HI 99-00 e+p incl. jets \59} 


1.5% 


0.9762 


0.9834 




ZEUS 96-97 e+p incl. jets \57} 


2% 


0.9735 


0.9811 




ZEUS 98-00 e±p incl. jets [58] 


2.5% 


0.9656 


0.9761 




D0 II pp incl. jets [56] 


6.1% 


0.9353 


1.0596 


1.0759 


CDF II pp incl. jets [54] 


5.8% 


0.8779 


0.9646 


0.9900 


CDF II ^ asym. [48] 




1 


1 


1 


00 111^^ £iy asym. [49] 




1 


1 


1 


D0 II Z rap. [53] 




1 


1 


1 


CDF II Z rap. [52] 


5.8% 


0.8779 


0.9646 


0.9900 



Table 3: The fitted normalisations J\fn of the data sets included in the global fit, together with 
the one-sigma normalisation errors, o"^, for each data set n. 
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where 



A 



kk' 



^kk' + 



pts. ^corr 



n,fc,i ^n,k',i 
^uncorr. 



i?A.({a},Ar„) = 5^ 



pts. corr. 



(D„,,-T„,,({a})/Ar„) 



uncorr. \ 2 
n,i J 



(40) 



Therefore, the optimal shifts of the data points by the systematic errors are solved for analyti- 
cally, while the input PDF parameters {a}, together with the data normalisations A/'„, must be 
determined by numerical minimisation of Xgiobai usual way. 

In practice, we use f p8|) only for the jet-energy-scale uncertainties in the ZEUS inclusive jet 



data |57l. l58l|. the 4 contributions to the correlated uncertainty in the HI inclusive jet data [59 



the 16 sources of correlated uncertainty in the CDF Run II inclusive jet data [5J], the 22 



sources of correlated uncertainty in the D0 Run II inclusive jet data [56|, and the 6 sources of 



correlated uncertainty in the CDF Run II Z rapidity distribution [52] ■ In many other cases the 



uncertainties (other than normalisation) are presented without any information on correlations 
between systematic uncertainties, or are assumed to be uncorrelated. For the remainder of the 
data sets, the systematic errors (other than the normalisation error) are simply treated as being 
uncorrelated and added in quadrature with the statistical uncertainties, i.e. we use with the 
penalty term (1571) . 



Minimisation of and calculation of Hessian matrix 

To determine the best fit at NLO and NNLO we need to minimise the Xgiobai with respect to 
28 free input PDF parameters, together with as{Qo), 3 parameters (ri, r2, r^) associated with 
nuclear corrections (see Section [7]), and 17 different data set normalisations, giving a total of 49 
free para mete rs. This is a difficult task and is unlikely to be possible with the widely-used minuit 
package 127l |. where the practical maximum number of free parameters is around 15. Indeed, 



the CTEQ group have found it necessary to extend minu it to use an improved iterative method 



in the numerical calculation of the Hessian matrix |128l |. Other fitting groups, such as Alekhin 
and Hl/ZEUS, use far fewer PDF parameters, and so presumably do not suffer from the same 
problems when using the standard minuit packag e. As in p revious MRS/MRST analyses |2|-|21 |. 
we instead use the Levenberg-Marquardt method 129l . ll3Cll | described in Numerical Recipes (l3l| , 



which combines the advantages of the inverse-Hessian method and the steepest descent method 
for minimisation. The method requires knowledge of the gradient and Hessian matrix of the 
Xgiobab which we provide partially as analytic expressions in our fitting code, only using numerical 
finite-difference computations for the derivatives of the theory predictions with respect to the 
fitted parameters. A linearisation approximation is made in the calculatio n of the Hessian matrix 
to avoid the presence of potentially destabilising second-derivative terms 131| . For example, the 
contribution to the global Hessian matrix Him 
(Bl is 



-^/m from a data set n with Xn defined by 



TTn — 
^Im — 



1 dxl 

2 dai dttr, 



pts. 



.a. 



n,i 



9T„,,({a})/Ar„ dT^,,{{a})IUn 



corr. \ 

n,kAJ 



dai 



(41) 
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while in the more comphcated case of treating correlated systematic errors (155]) , the correspond- 
ing expression is 

Tjn _ Sr^ 1 dTn^i dTn,i ST-^ dVn^k dfn,k , .n\ 

where T„_i({a}, AT^) = T„,i({a})/A/'n + Zlfcti" AT^X^^^. Similar expressions are used 

for the elements of the Hessian matrix corresponding to data set normalisations A/'„. 

5.3 Discussion of fit results 

For each of the three fits, the optimal values of the input PDF parameters of Eqs. fl6l)- f|T2l) and 
of the QCD couphng as are given in Table HI 

Notice that t he L O fit gives a significantly larger 05 (M|) than the world-average value of 



0.1176 ± 0.0020 |132l |. and that the parametrisation of the input gluon is simpler since the fit 
never shows any tendency for the gluon to go negative, or even to turn over at small x. Moreover, 
the LO fit gives a much worse description of the data than that obtained at NLO and NNLO. 
It fails for HERA structure functions because the evolution is too slow on account of missing 
higher-order small-x enhanced terms in the quark-gluon splitting function. It fails for some 
fixed-target data due to the absence of higher-order large-x enhanced terms in the coefficient 
functions. For the HERA jets it is missing large NLO corrections to the cross sections. For 
Tevatron jets we find that the LO evolution cannot give, in detail, the correct shape to fit the 
data, contrary to previous findings for the less precise Run I data which were fit very well It 
is also missing large NLO corrections to the cross section for the Drell-Yan data, and as before we 
have applied a i^-factor (previously fixed to 1.3), but have improved this to {l + as{M^) Cpn/2), 
which improves the fit quality considerably. The deficiencies due to higher-order terms cannot 
be completely mimicked by the very large coupling as (although this works quite well at high 
x) and a large small-x gluon, even when the normalisations of many data sets choose to lie 
rather low. So, in conclusion, the LO fit is unsatisfactory. On the other hand, if required, it is 
possible to obtain a modified LO PDF set which gives a much improved description of the data 
by relaxing the momentum conservation constraint and using the NLO coupling as |66j |. and 
this can be even further improved by a modification of the scale of the coupling [671]. 

The quality of the fit at NLO is much better, with no obviously badly fit data sets (except 
to some extent the D0 W —>■ iu charge asymmetry, which will be discussed later). Again the 
normalisations tend to lie a little low, though there are systematic differences between data 
sets, e.g. the SLAG DIS data are normalised about 1.5% higher than NMC, which is 3% higher 
than BCDMS, and in general HI sets are 1% or so higher in normalisation than ZEUS. These 
trends have long been known. We also see that the fitted normalisation of the D0 Run II jet 
data is about 8% higher than that for CDF Run II data (in the latter case the normalisation of 
the inclusive jet data is tied to the normalisation of the Z rapidity distribution data), but the 
luminosity uncertainties in each are compatible with this. The fit quality at NNLO is almost 
identical to that at NLO, and is rather similar for most data sets. However, it is striking that 
the normalisations are clustered about unity much more symmetrically in this case, with those 
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Parameter 


LO 


NLO 


NNLO 


as{Ml) 


0.68183 
0.13939 


0.49128 
0.12018 


0.45077 
0.11707 


A-u 

m 

7« 


1.4335 
0.45232 

3.0409 t°:r7 

-2.3737 t^^f^ 
8.9924 


0.25871 
0.29065 

q 9/1 qo +0.062 
o.zioz. -0.039 

4.0603 
30.687 


0.22250 
0.27871 

q qG07 +0.061 
O.ODZ/ _o.044 

4.4343 ^2:7 
38.599 


Ad 

V3 

m - m 

Id 


5.0903 
0.71978 

2.0835 
-4.3654 

7.4730 


12.288 
0.96809 t^ii 
2.7003 +°-^° 

-3.8911 

6.0542 


17.938 

1.0839 +1^11 
2.7865 
-3.6387 
5.2577 


As 
Ss 
Vs 
es 

7s 


0.59964 
-0.16276 

8.8801 t":! 
-2.9012 
16.865 


0.31620 t'^;]^^ 
-0.21515 

9.2726 t'^^-'l 
-2.6022 
30.785 


0.64942 

-0.11912 

9.4189 +^3^ 
-2.6287 
18.065 


/;dxA(x,Q§) 
Aa 

VA 
7A 

Sa 


0.091031 
8.9413 
1.8760 

8.4703 tg!^ 
-36.507 


0.087673 tg-^i? 

8.1084 

1.8691 

13.609 
-59.289 


0.078167 
16.244 

2.0741 +0:18 
6.7640 tVl 
-36.090 


Ag 

s, 

Vg 
^9 

7s 

Ag' 

Sg' 

Vg' 


0.0012216 

-0.83657 

l.Ami „o.5o 

-38.997 
1445.5 


1.0805 
-0.42848 1°;°^? 
3.0225 -0.36 
-2.2922 
3.4894 
-1.1168 
-0.42776 
32.869 


3.4055 
-0.12178 

2.9278 to;4i 
-2.3210 
1.9233 
-1.6189 
-0.23999 ^;^;io 
24.792 2 


A+ 

V+ 
A— 

V- 

Xo 


0.10302 +^:;jf^ 

13.242 +1-^ 

n ni 1 t;oQ +o!oo9 
-0.011523 _ooi8 

10.285 

0.017414 


0.047915 

9.7466 tj-g 

n ni 1 Ron +0.009 
-0.011629 _oo23 

11.261 

0.016050 


0.10455 

9.8689 tJ:6 

n nnno^^no +0.006 
— U. 0093692 1q 024 

9.5783 +f 
0.018556 


r\ 

Tl 

r-i 


-0.39484 
-1.0719 
-0.28973 


-0.57631 
0.81878 
-0.083208 


-0.80834 
1.2669 
0.15098 



Table 4: The optimal values of and the input PDF parameters at = 1 GeV^ determined 
from the global analysis. The one-sigma errors are calculated using flSip and (1521) using the 68% 
C.L. tolerance discussed in Sectional and are shown only for the 20 parameters allowed to go free 
when determining the eigenvector PDF sets. The parameters A^^ Ad, Ag and xq are determined 
from sum rules and are not fitted parameters. Similarly, A^ is determined from dx A{x, Qq)- 
The three parameters r^, defined in (1731) . are associated with the nuclear corrections to the 
neutrino data; see Section I7.3[ The parameter values are given to five significant figures solely 
for accuracy in the case of reproduction of the PDFs. 
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for the E866/NuSea Drell-Yan cross sections and the D0 Run II jet data more than one-sigma 
above unity. At NLO and NNLO the values of aJMl) of 0.1202 and 0.1171 are both consistent 
with the world average value of 0.1176 dz 0.0020 132| . though the NLO value is just a little high, 



suggesting a partial compensation for missing higher orders. A detailed study of the uncertainty 
in as awaits a further publication, but previous estimates of an uncertainty on our determination 
of asiM"^) of ±0.002 from data and ±0.003 from theory are unlikely to change very significantly. 

The errors that are shown on the input parameters in Table H] are calculated with (15T|) 
and (152|) using the 68% C.L. tolerance discussed in Section [6], and are given only for the 20 
parameters allowed to go free when determining the eigenvector PDF sets. It is important 
to note that although they give some idea of PDF uncertainties they are only quantitatively 
meaningful when the information on the correlations between parameters is included, and this 
is most easily presented and utilised in terms of orthogonal eigenvectors. The eigenvector PDF 
sets are defined and discussed in Section 16. 1[ Indeed, in Section [6] we describe, in detail, the 
procedure that we use to determine the uncertainty on the parton distributions and on the 
predictions of cross sections. We then turn to the more detailed description of the various data 
sets that are obtained in the global fits in Sections [71-[T2l 
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6 Error propagation from experimental data points 



Uncertainties in global PDF analyses can be divided into two general categories, which may 
loosely be called "theoretical" and "experimental" . The theoretical errors include those associ- 
ated with the choice of the form of the input parameterisation, the neglected higher-order QCD 
(including enhanced ln(l/a;) or ln(l — x) terms) and electroweak corrections, parton recombina- 
tion and other higher-twist corrections, the choice of data sets and kinematic cuts, the choice of 
nuclear corrections for the neutrino-initiated data and the treatment of heavy flavours. These 
uncertainties are often difficult to quantify a priori until a better calculation or prescription 
becomes available, but see the attempt made in Ref. [ITJ. 

On the other hand, in principle there are well-defined procedures for propagating experimental 
uncertainties on the fitted data points through to the PDF uncertainties. Three main methods 
have been used to do this. 



1. The Lagrange multiplier method |128l . ll33l |. which does not rely on linear error propagation, 
but requires the ability to perform a global fit. 



2. The Hessian method |134| . which is based on linear error propagation and involves the 
production of eigenvector PDF sets suitable for convenient use by the end user. 



3. The use of Monte Carlo sampling [l35l . Il36| . which has recently be en us ed in conjunction 
with neural networks to determine NLO PDFs from a DIS-only fit 137 . 



The first two methods were originally used by CTEQ |128l . Il33l . Il34 | and then by MRST [16 



In principle, the Lagrange multiplier method is superior to the Hessian approach, but it suffers 
from the enormous practical disadvantage that a series of new global fits has to be done every 
time one considers a new quantity. Fortunately, it turns out that for those quantities that have 
been considered by both methods, the uncertainties have been found to be comparable. We 
therefore concentrate here on the most commonly used Hessian method. We will not discuss the 
use of Monte Carlo sampling in this paper. 



6.1 Review of Hessian method in global PDF analyses 

The basic procedure is discussed in detail in Refs. fiel-fii^. Here, we briefiy review the impor- 
tant points. The traditional propagation of experimental uncertainties assumes that the global 
goodness-of-fit quantity, Xgiobai; quadratic about the global minimum, which has parameters 
{a?, . . . , a°}. In this case we can write 

n 

^Xgiobai = Xgiobai - Xmin = ^ Hij{ai - a°)(aj - a°), (43) 



where the Hessian matrix H has components 



-L A global 



2 doidoj 



(44) 
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The uncertainty on a quantity F{{ai}) is then obtained from hnear error propagation: 



AF 



—c — 

^ dai da A 

i,j=i ■> 



(45) 



where C = is the covariance matrix, also known as the error matrix, and T — v^Agiobaiy' 
is the tolerance for the required confidence interval. This formula fHSj) has the disadvantage that 
PDF uncertainties are not readily calculable for general observables, since the derivative of the 
observable F with respect to each parameter is needed. 



It is convenient to diagonalise the covariance (or Hessian) matrix 16l . Il34| . and work in terms 
of the eigenvectors and eigenvalues. Since the covariance matrix is symmetric it has a set of 
orthonormal eigenvectors Vk defined by 



CijVjk 



AfcW 



ik 7 



(46) 



where is the kth eigenvalue and Vik is the ith component of the kth orthonormal eigenvector 
{k = 1, . . . ,n). The parameter displacements from the global minimum can be expanded in a 
basis of rescaled eigenvectors 6,^ = v^Afefjfc, that is. 



tti - a, 



E 

k=l 



^ik^k- 



(47) 



Then it can be shown, using the orthonormality of Vk-, that (H3l) reduces to 



A global 



Amin 



k=l 



■ki 



(4J 



that is, 'Y^=i ^k — is interior of a hypersphere of radius T. Pairs of eigenvector PDF 
sets can then be produced to span this hypersphere, at a fixed value of as, with parameters 
given by 

^±^-a°±te,fc, (49) 

)^/^. In the quadratic approximation, t = T. 
For the larger eigenvalues A^, where there are significant deviations from the ideal quadratic 
behaviour, t is adjusted iteratively to give the desired value of T. Then uncertainties on a 
quantity F, which may be an individual PDF at particular values of x and Q^, or a derived 
quantity such as a cross section, can be calculated with0 



with t adjusted to give the desired T = (AXgiobaU 



AF 



1 



(50) 



k=l 



^^It can be shown that reduces to (j5D)) in the quadratic approximation {t — T) [l^, 134 1. but we treat 
(|50p as the fundamental definition in the departure of this ideal limit. In this paper, we will generally use (|5ip 
and (|52p to calculate asymmetric PDF uncertainties. 
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or asymmetric errors can be calculated with 



(AF)_ 



J2 {max [ - F(5o), F{S,) - F(5o), O] }' 

\ k=l 



J2 {max [ FiSo) - F(5o) - FiS,), O] }^ 

\ fc=i 



(51) 



(52) 



where 5*0 is the central PDF set. Correlations between two quantities can also be calculated; 
see, for example, Ref. 87!]. Defining a correlation cosine between two quantities F and G, 



cos (pFG 



4AFAG 



k )\ ' 



(53) 



k=l 



where the uncertainties AF and AG are calculated using fl5U]) . then values of cos(f)FG ~ 1 mean 
that F and G are correlated, values of ~ — 1 mean that they are anticorrelated, while values of 
« mean that they are uncorrelated. A tolerance ellipse in the F-G plane can then be defined 
by the two parametric equations: 



F = F{So)+AF cose, 

G = G{So) + AG cos(^ + (f)FG) 



(54) 
(55) 



where 9 e [0,27r]. We will show examples in Section [T5l Note that u sing e.g. 90% C.L. PDFs 



results in a probability less than 90% for the 2-D tolerance ellipse 87|, Il3ll . Il32 . 

To determine the "best-fit" parameters we allow all the input PDF parameters of Eqs. 



(jT2|) to vary. However, when investigating in detail the small departures from the global minimum 
we notice a certain amount of redundancy in parameters. A striking example is for the NLO 
and NNLO parameterisations of the gluon distribution flTU]) . Small changes in the values of 
three of the parameters can be compensated almost exactly by changes in the remaining four. 
This high degree of correlation between parameters means that very small changes in will 
be obtained. However, at some point the compensation starts to fail significantly and the 
increases dramatically. Hence, the redundancy leads to a severe breaking of the quadratic 
behaviour in Ax^, and some very fiat directions in the eigenvector space (that is, very large 
eigenvalues of the covariance matrix) and cubic, quartic etc. terms dominate. During the process 
of diagonalisation this bad behaviour feeds through into the whole set of eigenvectors to some 
extent. Therefore, in order that the Hessian method works at all well we have to lessen the 
redundancy in the input parameters. In order to do this we simply fix some of the parameters 
at their best-fit values, so that the Hessian matrix only depends on a subset of parameters that 
are sufficiently independent that the quadratic approximation is reasonable. We finish up with 
20 reasonably well-behaved eigenvectors in total, i.e. those coming from the combinations of 
the 20 parameters that are assigned errors in Table |H However, we emphasise that the other 
parameters are fixed at their best-fit values, rather than simply set to zero. The problem of a 
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Figure 5: The points (•) show A^gj^j^^^j as a function of the distance along each eigenvector 
direction, t, defined in ( l49l) . for eigenvectors numbered 1-10 corresponding to the 10 smallest 
eigenvalues. The dashed curve is the ideal case, AXgi^^j^^^j = t^. 
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Figure 6: The points (•) show A^gj^j^^^j as a function of the distance along each eigenvector 
direction, t, defined in ( l49l) . for eigenvectors numbered 11-20 corresponding to the 10 largest 
eigenvalues. The dashed curve is the ideal case, Ax^i^^j^j^j = t^. 
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certain amount of redundancy in parameters is a general feature of the full global fits obtained 
by CTEQ and MRST which have sufficient parametric freedom to account for all features of the 
data. 

The behaviour of Ax^j^j^j^j for the 20 eigenvectors in the current NLO fit is shown in Figs. [5] 
and El For the lowest 13 or so eigenvalues the quadratic approximation is very good. Above 
this we see some departure from quadratic behaviour, though there is still good agreement with 
the ideal case for relatively low t. The peculiar behaviour for eigenvector 20 is related to the 
fact that this is mainly concerned with the high-x part of = s — s at input, and this will be 
discussed in Section 16.51 



6.2 Dynamic determination of tolerance 

Ideally, with the standard "parameter-fitting" criterion, we would expect the errors to be given 
by the choice of tolerance T = (Ax^i^^j^j^j)^/^ = 1 for the 68% (one-sigma) confidence level (C.L.) 



limi 



or 



LglobaL 

2.71 for the 90% C.L. limit 



131, 132 



This is appropriate if fitting consistent 
data sets with ideal Gaussian errors to a well-defined theory. However, in practice, there are 
some inconsistencies between the independent fitted data sets, and unknown experimental and 
theoretical uncertainties, so the parameter-fitting criterion is not appropriate for global PDF 
analyses. One possible app roach is to attempt to only include completely consistent data sets. 



This was tried in Ref. |136l |. but it was found that the ZEUS and NMC data were inconsistent 
with HI, BCDMS and E665, so only the latter three data sets were used in their fit. Clearly, 
this does not seem like a practical approach to follow. 

So how can we choose the value of the tolerance T so as to obtain more reliable errors? 
Instead, we can appeal to the much weaker "hypothesis-testing" criterion, where the eigenvector 
PDF sets are treated as alternative hypotheses. Very roughly, a fit is judged to be "good" if 



138(1 . More precisely, ranges 



each data set n, consisting of data points, has Xn — =t V'2N 
of Xn corresponding to a 90% C.L. limit, for example, can be calculated, then the value of the 
tolerance T = (AXgi^^^^j^j)^/^ can be chosen to ensure that each data set is described within its 
90% C.L. limit. Assuming that Xn follows the x^-distribution with A^ degrees of freedom, which 
has a probability density function 



Pn{x' 



(x 



2^Af/2-l g-xV2 



2JV/2r(]v/2) 

then the 90th percentile, ^90, is obtained by solving 

rS.90 



(56) 



/ dx' Pwix') = 0.90. 
Jo 



(57) 



Similarly for the 68th percentile, ^qs, and for the most probable value, ^50 — A^. In Fig. [7] we 
show (156!) for A^ = 163 corresponding to a typical data set included in the global fit (BCDMS 
fxp F2). We have ^so = 162.3, ^68 - ^50 = 8.7 and ^90 - ^so = 24.2, cf. V2N = 18.1. 



^^To be precise, rather than using exactly 68% we use er f(l/V2) = 0.682689. 
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probability density function for N = 163 d.o.f. 




Figure 7: Probability density function of the x^-distribution with 163 degrees of freedom. 

We now describe the new procedure we adopt in determining the tolerance for each eigenvec- 
tor direction: it f ollows clo s ely, but significantly extends, similar studies carried out earlier by 



the CTEQ group |l26l . Il33l . Il34j . Since the new procedure is quite complicated, we will break 



the explanation down into steps, supplemented by example plots from the current NLO analysis. 

(i) First we plot the change in Xn each data set, as the difference from the value at the 
global minimum, Xnc when moving along a particular eigenvector direction in units of 
i^XgiohaiY^'^- example for a subset of the most constraining data sets for eigenvector 
number 13 is shown in Fig. [HI The points (•) are generated for fixed values of the distance, 
i^XgiohaiY'''^^ between and 10 in each "+" or "— " eigenvector direction. Note that we 
adjust t in fH9l) until the desired i^XgiohaiY^'^ obtained: these two quantities are equal 
only in the strict quadratic approximation. These points (•) are then fitted to a quadratic^ 
function shown by the solid lines. Note from Fig. [S] that the minimum value of Xn for a 
particular data set does not in general occur at the global minimum, indicating some 
tension between data sets included in the global fit, for example, between the E866/NuSea 
pp Drell-Yan data and the NuTeV and CHORUS uN xF^ data. 

(ii) We define the 90% C.L. region for each data set by the condition that 

Xl < (1^) ^90, (58) 
and similarly for the 68% C.L. region. Note t hat we define the 90% C.L. region (l58l) for 



by renormalising ,^90 by a factor Xno/^so [l33|. This procedure is necessary to take 



""^•^The one exception is for the very asymmetric x\ profile for the NuTeV dimuon data plotted against the 
distance along eigenvector number 20, associated with the 77^ parameter controlling the high-x behaviour of 
Sv = s — s, which we instead fit to a quartic function. 
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Figure 8: Change in Xn for a subset of data sets n, compared to the values at the global minimum 



2 

n,0 



), when moving along eigenvector number 13 in units of {'^^Xgiohai 
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C.L. regions determined according to fl58l) are indicated. 
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account of the fact that the value of Xn o global minimum may be quite far from the 

most probable value of ^50 of the particular data set n. For the example of the BCDMS 
HP F2 data, the value of the Xn our NLO fit is 182.2, which would lie outside the 68% 
C.L. region shown in Fig. [7] if the absolute values of Xn were used instead. After applying 
the rescaling factor of Xno/^so = 1-12, the 90% C.L. region is given by Xn ~ Xno < ^7.2 
and the 68% C.L. region hj Xn~ xio < The 90% and 68% C.L. regions determined 
in this way are shown as the horizontal dashed lines in Fig. [HI 

(iii) For each eigenvector, we then record the value of the distance, (AXgiobai)^^^' which the 
Xn for each data set is minimised, together with the 90% and 68% C.L. limits defined by the 
intercepts of the quadratic curves with the horizontal dashed lines in Fig. [HI These values 
are shown for eigenvectors 9 and 13 in Fig. [9], where the points (•) indicate the values of 
the distance, (AXgiobai)^^^' fo^ which Xn is minimised, while the inner error bars extend 
across the 68% C.L. region and the outer error bars extend across the 90% C.L. region 
defined by 



(iv) For each of the eigenvectors, we choose the values of the tolerance T = (Axgj^j^^g^j)^/^, 
indicated by the horizontal dashed lines in Fig. [HI so that all data sets are described within 
their 90% or 68% C.L. regions. For eigenvector 9 we see that the tolerance in the "+" 
direction is fixed by the HI NC data and in the "— " direction by the ZEUS NC data, 
while for eigenvector 13 we see that the tolerance in the "+" direction is fixed by the 
E866/NuSea pp Drell-Yan data and in the "— " direction by the NuTeV uN xF^ data. 

We summarise the values of the tolerance obtained by this procedure in Fig. [10] for each of the 
20 eigenvectors of the NLO fit. The inner error bars indicate the tolerance for a 68% C.L. limit 
while the outer error bars indicate the tolerance for a 90% C.L. limit. The labels placed at the 
end of the error bars indicate the name of the data set which fixes the value of the tolerance 
for a 90% C.L. limit. Note that the results shown for eigenvectors 9 and 13 correspond to the 
values that we extracted from the plots of Fig. [9] 

The procedure just outlined is an extension of that previously employed by CTEQ and 
MRST. These previous analyses used a fixed value of the tolerance T = {^X^^ioh&x)^^'^ producing 
eigenvector PDF sets, namely T = vTOO for CTEQ and T = v^50 for MRST. Analogous plot s 



to those shown in Fig. [9] were originally presented by CTEQ in Appendix B.4 of Ref. [126 . 
Inspection of Fig. [TOl shows that our new dynamic tolerance required to ensure that all data sets 
are described within their 90% C.L. limits are almost all in the region T ~ a/50, which is close 
to the MRST value and suggests that the CTEQ tolerance (T = vToO) is too large0 However, 
an even smaller tolerance is obtained for some eigenvectors, in particular those associated with 
strange quarks and the d — u difference, both of which are constrained by a small number of 
data sets generally each with a relatively small number of data points. 



^ Note, however, that recent CTEQ analyses [8J,[87| apply unspecified weight factors to the data from certain 
experiments in their Xgiobai definition, which in practice achieves a similar effect to our dynamic tolerance method. 
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Figure 9: Ranges of (AXgiobai)^''^ along (a) eigenvector 9 and (b) eigenvector 13 for which data 
sets are satisfied within their 90% C.L. hmit (outer error bars) or 68% C.L. hmit (inner error 
bars). The points (•) indicate the minimum with respect to each particular data set. The 
tolerance, indicated by the horizontal dashed lines, is chosen to ensure that all data sets are 
described within their 68% or 90% C.L. limits defined by 
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MSTW 2008 NLO PDF fit 
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Figure 10: Tolerance for each eigenvector direction determined dynamically from the criteria 
that each data set must be described within its 90% C.L. (Eq. (158|) ) (outer error bars) or 68% 
C.L. limit (inner error bars). The labels give the name of the data set which sets the 90% 
C.L. tolerance for each eigenvector direction. 



6.3 Uncertainties on input PDFs 

We use the values of the dynamic tolerance shown in Fig. [10] to generate the PDF eigenvector 
sets according to (H9l) . which can be written as 

a,(5±) = a°±t±e,,, (59) 

where is adjusted to give T^, with the values shown in Fig. [101 We provide two different 
sets for each fit corresponding to either a 90% or 68% C.L. limit. Note that the ratio of the PDF 
uncertainties calculated using these two sets is not simply an overall factor of y/2.71 = 1.64, as 
it would be if choosing the tolerance according to the usual parameter-fitting criterion. Even 
in the simplest case, where the data set fixing the tolerance is the same for the 90% and 68% 
C.L. limits, and assuming linear error propagation, then the ratio of the values would be 
(^90 ~ C5o)/(^68 ~ ^5o)) which is a function of the number of data points in the data set which 
fixes the tolerance, and takes a value around 1.7 for typical ~ 10-1000. 
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MSTW 2008 NLO PDF fit (68% C.L.) 

Contribution to uncertainty from eigenvector number 9 
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X 

Figure 11: The total uncertainty on the input PDFs (shaded bands), given by flSTj) and fl52l) . 
and the contribution from = 9 in these equations (dashed hues). For this eigenvector, the 
dashed hues are hidden under the central solid lines with the exception of xg and xs + xs. 
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The total one-sigma uncertainty on the input parton distributions at Qq = 1 GeV^, calcu- 
lated using fl5T|) and fl52|) . is shown as the shaded bands in Fig. [TTl Note that the somewhat 
peculiar shape of the uncertainty on s — s is clearly due to the limited freedom in the input pa- 
rameterisation: only two free parameters controlling the normalisation and the power of (1 — x). 

Recall that each eigenvector Vk is a linear combination of the input PDF parameter dis- 
placements tti — a° and that the corresponding eigenvector PDF sets 5*^ are responsible for a 
certain uncertainty on the parton distributions. Sometimes an eigenvector corresponds to one 
particular input distribution, but often a single eigenvector contributes to the uncertainty on 
a variety of the input distributions. The dashed lines in Fig. [TT] show the contribution to the 
uncertainty from a single eigenvector (in this case, number 9), i.e. only the k = 9 term in ( ISTl) 
and (l52l) . Clearly, eigenvector 9 mainly contributes to the uncertainty on the input gluon distri- 
bution. This can be made clearer by plotting the fractional contribution to the uncertainty on 
the input PDFs from a single eigenvector, i.e. the contributions to fl5T|) and fl52|) from a single 
value of k divided by the total uncertainties. The fractional contribution to the uncertainty from 
eigenvector number 9 is shown in Fig. [12] and that from eigenvector number 13 in Fig. [131 

A general summary of all the plots of the type shown in Figs. [12] and [13] is presented in Table 
[5] The bold entries show which eigenvectors k are mainly responsible for the uncertainty on 
the various input distributions. The three numbers in each entry correspond to (ten times) the 
fractional contribution to the uncertainties (AF^/AF for the F = g,Uv, . . . input distributions, 
where AFk is the contribution to ( ISTl) and ( |52l) from a single eigenvector k) in the small x 
{x < 0.01), medium x (0.01 < x < 0.1) and large x {x > 0.1) regions, respectively, for the NLO 
global fit with 68% CL. uncertainties. We see that each input distribution is sensitive to only 
a small number of eigenvectors. An indication of which data sets are particularly important for 
which flavour of parton distribution can be obtained by correlating the results of Fig. [10] with 
the table. For example, from Table [5] we see that the gluon distribution is particularly sensitive 
to eigenvectors 9, 11 and 19. The tolerance for each of these eigenvectors is fixed by the HI and 
ZEUS NC data, apart from for the 68% CL. tolerance for eigenvector 19 in the "— " direction, 
where it is fixed by the CDF Run II high-pj- jet data. 

In fact, the entries in Table [5] can be used in conjunction with the eigenvector plots, such as 
those used as examples in Fig. [9] to pinpoint which particular data sets constrain the various 
parton distributions. For example, the entries in Table [5] show that eigenvector 13 gives the 
biggest contribution to the high-x sea distribution uncertainty and is also very significant for the 
d — u distribution at high x. Inspection of the plot of eigenvector 13 versus the data. Fig. Mjo), 
shows that the tightest constraints on the eigenvector 13 come from the data sets collected in the 
following experiments: E866/NuSea pp and pd/pp Drell-Yan, NuTeV and CHORUS uN F2 and 
NuTeV and CHORUS uN xF^. The former are the obvious direct constraint on the sea quarks, 
while the others influence the decomposition into the well-determined valence quarks and sea 
quarks. Eigenvector 13 also contains information on the dy uncertainty, which is determined 
by the neutrino DIS data which weight down quarks more highly than neutral-current DIS on 
proton targets. Similarly, eigenvector 9 only influences the gluon uncertainty to any great degree, 
so it is no surprise that it is constrained by the HI and ZEUS NC data (and CDF Run II high-pj- 
jet data); see Fig. [Hl^a). Surprisingly, the summary plot in Fig. [TU] shows that the Tevatron jet 
data are never the main constraint on any eigenvector, apart from for the 68% CL. tolerance for 
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MSTW 2008 NLO PDF fit (68% C.L.) 

Fractional contribution to uncertainty from eigenvector number 9 










At input scale 
Q2 = 1 GeV2 



Figure 12: The fractional contribution to the uncertainty on the input PDFs from eigenvector 
9, defined by the contributions to (ISTIl and (l52l) from k = 9, divided by the total uncertainties. 
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MSTW 2008 NLO PDF fit (68% C.L.) 

Fractional contribution to uncertainty from eigenvector number 13 










At input scale 
Q2 = 1 GeV2 



Figure 13: The fractional contribution to the uncertainty on the input PDFs from to eigenvector 
13, defined by the contributions to (!5T]) and (l52l) from k = 13, divided by the total uncertainties. 
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Table 5: The three numbers in each entry are the fractional contributions to the uncertainties 
{AFk/AF for the F = g,Uy, . . . input distributions, where AF^. is the contribution to fl5T]) and 
( 152|) from a single eigenvector k) in the small x {x < 0.01), medium x (0.01 < x < 0.1) and 
large x {x > 0.1) regions, respectively, arising from eigenvector k in the NLO global fit, using 
68% CL. errors. Each number has been multiplied by ten; for example, 4 denotes 0.4. Major 
contributions are shown in bold type. For a precise value of x, and symmetric errors, the sum 
of the squares of each column should be 100. However, the entries shown are the maximum 
fraction in each interval of x, so often do not satisfy this condition. 



eigenvector 19 in the "— " direction (not shown here). This does not necessarily mean that the 
Tevatron jet data provide no constraint, however, and this will be discussed in detail in Section 

A complete set of plots related to the error analysis and PDF uncertainties discussed in 



Sections 16.21 and 16.31 similar to the examples shown in Figs. [5]-[T3l is available at Ref. [27 . 
These supplementary plots should prove valuable for the end users of the PDFs, in investigating 
the sensitivity of different observables to different eigenvector directions, and in finding the 
corresponding data sets which constrain those particular eigenvector directions. 
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6.4 Effect of free data set normalisations on PDF uncertainties 



As already noted, one of our significant improvements in the uncertainty determination in tfie 
present analysis is the inclusion of the data normalisation errors while determining the PDF 
uncertainties. In previous MRST (and CTEQ) fits these have been held fixed during the error 
propagation. To account for free data set normalisations, we first calculate the 39 x 39 covariance 
matrix with an additional 19 free parameter^ corresponding to data set normalisations, in 
addition to the usual 20 free input PDF parameters. In principle the inclusion of free data set 
normalisations in PDF uncertainties can be done in two different ways. 

The simplest way is to diagonalise the 39 x 39 covariance matrix to obtain 39 eigenvectors 
and to determine the tolerance values dynamically as described above. However, this method 
involves a large number of eigenvectors, many of which actually consist largely of the variation of 
data set normalisations. It is clearly more physically intuitive to use eigenvectors corresponding 
only to input PDF parameters. 

An alternative method, which has this advantage and which we use in the present analy- 
sis, is to remove the rows and columns from the 39 x 39 covariance matrix corresponding to 
data set normalisations, leaving only a 20 x 20 covariance matrix corresponding to input PDF 
parameters, which contains, however, information on free data set normalisations from the in- 



version of the Hessian matrix [131| . Then the 20 x 20 covariance matrix can be diagonalised to 
give 20 eigenvectors, and the tolerance values and eigenvector PDF sets produced as described 
above, with the modification that every time the PDF parameters are changed, when moving 
along an eigenvector direction, the Xgiobai should be minimised with respect to the data set 
normalisations!^ 

In the limit of quadratic behaviour for the Xgiobai function and equal tolerance for each 
eigenvector we would obtain the same results for PDF uncertainties in each case. In Fig. [H] we 
show the effect of allowing data set normalisations to go free according to these two methods, 
compared to the case where data set normalisations are fixed as in the MRST analyses. The 
new procedure leads to an increase in uncertainties, but in most places this increase is very 
small. The exception is those parton distributions where the uncertainty is smallest, e.g. the up 
valence quark for x > 0.05 and the light sea for x between 0.001 and 0.01. The absolute increase 
in uncertainty is 1% or less, but in these cases it is a significant contribution to the total. We 
see that the approach using data set normalisations in eigenvectors gives very similar results, 
demonstrating that despite our quartic penalty term for normalisations there is little breakdown 
of the quadratic approximation. Indeed, the uncertainty on normalisations determined from the 
fit and applying the tolerance approach is usually rather similar to the quoted experimental 
uncertainty. Hence, we do not actually encounter large penalties from varying normalisations, 
as we rarely vary by much more than the experimental one-sigma uncertainty. 



"'^^To be precise, there are 16 free parameters associated with data set normalisations together with the 3 
parameters (ri, r2, r^) associated with nuclear corrections ((73|) . We do not allow the normalisation of the neutrino 
data to go free because it is very highly correlated with ri . 

^^Similarly, the Xgiobai should be minimised with respect to the shifts in correlated systematic errors, in 
Eq. each time the PDF parameters are changed, but this is done automatically by using Eq. 
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Figure 14: The fractional uncertainty of the MSTW 2008 NLO PDFs at = QeV^ us- 
ing 20 eigenvectors and allowing for free data set normalisations. We also show the fractional 
uncertainty obtained using 20 eigenvectors with fixed data set normalisations (as in the MRST 
analyses) and with 39 eigenvectors including 19 parameters associated with data set normalisa- 
tions in additional to the 20 input PDF parameters. 
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6.5 Input parameterisation and uncertainties 

There is a strong relationship between the input PDF parameterisation and the uncertainties 
which will be obtained. Too restrictive a parameterisation can limit the calculated uncertainty 
in an artificial manner. The parameterisation for the input PDFs in our analysis were presented 
in Section [31 specifically in Eqs. (l6i)- (fT2l) . The 28 free PDF parameters listed there allow a large 
degree of fiexibility, but are sufficiently few that only a few turning points are possible for each 
distribution, presumably the case in reality for quantities closely related to physical variables. 
Indeed, as outlined in Section 16.11 many of the 28 free PDF parameters turn out to be very 
strongly (anti) correlated, and only 20 parameters remain after removing all combinations where 
this is the case, whose uncertainties are quoted in Table HI This is 5 more parameters than in 
the MRST error sets |21| . Here we briefiy discuss the truly free parameters for each parton 



flavour, i.e. those parameters allowed to go free when calculating the covariance matrix used for 
error propagation, and outline the resulting effect on the uncertainty. 

Valence quarks 

For the up and down valence quark distributions we use the standard "MRS-type" parameteri- 
sation, 

xv{x, Ql) = A x"^" (1 - x)''" (1 + e„ + 7vx), (60) 

where the normalisation parameter is determined from the number sum rules (|T^ . In both 
cases the 3 parameters left free in the determination of the eigenvectors are 6y, rj^ and e„. These 
correspond to three reasonably independent regions of x. At high x the uncertainty is determined 
very largely by rj^. In detail, 

v±Av ~ A (1 + + 7,) (1 - x)""^^'''' 

~ vil-x)^^"!^ (61) 
~ V [iT^Vv ln(l - x)], 

where the last line applies if {Arj^ ln(l — x)\ is small enough for the linear expansion of the 
exponential to be valid. Hence, the parameterisation allows the uncertainty to become very 
large required by the lack of data in this regio n. I t is not as flexible as some 



of the competing parameterisations, e.g. the CTEQ6M PDFs |126l |. However, it does limit the 
possibility of unusual shapes at very high x, or arguably undesirable features, e.g. the gluon 
becoming the hardest distribution at very high x {x > 0.9), as found in the CTEQ6.1M set. At 
small X the uncertainty is determined very largely by 6^. Again in more detail, 

x{v±Av) ~ AyX^^^^^^ 

~ XV x^^^^ (62) 
~ xv[l±A6y\n{l/x)], 

where again the last line applies if |A(5t, ln(l/x)| is small enough for the linear expansion of the 
exponential to be valid. Similarly to the case at high x, the uncertainty can become very large 
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Figure 15: The fractional uncertainty on the (a) up valence quark and (b) down valence quark 
distributions at = iqO GeV^ from the MSTW 2008 NLO fit compared to that from the 
MRST 2001 NLO fit M 



at small x, but once the behaviour has been determined by the data at fairly low x, fixing the 
one parameter, even lower x values are simply an extrapolation. The third free parameter, e^, 
controls the variation of the input distribution at intermediate values of x {jv is more strongly 
correlated to rj^). In practice the uncertainties at very high x, and even more so at very small 
X, are also constrained by the sum rules. In the region where valence PDFs are constrained by 
data, i.e. x = 0.01-0.75, about 75% of the total number of quarks are found, so the number 
found in the extrapolation regions is relatively small and very well determined. 

The fractional uncertainty for the up valence distribution at NLO is shown in Fig. [T5l(a) 



along with that from the MRST 2001 analysis [16|. In this case the parameters left free when 



determining the eigenvectors are the same and the uncertainties are larger mainly because of the 
fact that data set normalisations were fixed in the MRST analysis; see Fig. HH The fractional 
uncertainty on increases rapidly for very high x and also for x < 0.003. In the latter case 
the data constraint for x < 0.01 at lower has moved to a smaller x due to evolution. The 
uncertainty for the down valence distribution is shown in Fig. [T5l(b). The free parameters in 
the MRST 2001 analysis were rja, and 7^. Hence, the MRST 2001 input parameterisation 
allowed less flexibility in the uncertainty at small x and from the sum rule this fed in to the 
unconstrained high-x region also. This is clearly illustrated in Fig. [TMb). 

Gluon 

The gluon distribution is the most complicated case for PDF uncertainties and parameterisations. 



The input gluon in the original MRST fits [12, ll3j had the same form of parameterisation as 



Eq. (!60|) . with Ag determined by the momentum sum rule. However, with the advent of improved 
HERA data included in the MRST 2001 fit [l^, it was noticed that if the evolution is started 
from Qg = 1 GeV^ the input gluon preferred to be negative at very small x and an additional 



56 



term Agi x^a' (1 — x)^^' was added to facilitat e thi s. To this date all competing parameterisations 



apart from the very recent NNPDFl.O set |137l) I . have a gluon which at small x behaves like 



xg{x, Ql) ~ x^^ , i.e. is controlled by a single power. This means that 

g±^g^g[l±MM^I^)l (63) 

i.e. the uncertainty grows linearly with ln(l/x) and there is no scope for a rapidly expanding 
uncertainty as data constraints run out. This is much more of an issue for the gluon than for 
valence quarks, as the momentum sum rule offers a far less direct constraint than the number 
sum rules as a; — >^ 0. However, there is another complication to consider, namely 

A^7(x, Ql) ~ Ql) A6g ln(l/x), (64) 

and so as g{x,Ql) becomes smaller then so does Ag{x,Ql). If g{x,Ql) is very small, then 
the absolute input uncertainty for the gluon is very small, and at higher the uncertainty 
is therefore determined entirely by evolution from higher-a;, i.e. by the region where the gluon 
distribution is better determined. Most PDF fitting groups find that xg{x, Q^) is indeed small 
at low and small x. In this region the MRST (since 2001) and MSTW gluon distributions 
have the form, 

xg{x,Ql) = xgi{x,Ql) +xg2{x,Ql) ~ AgX^' + Ag, x^a' ^ (65) 

which is more flexible than a single power. Not only does it allow the gluon to become negative 
at very small x, but it is also particularly important for the uncertainty, 

Ag{x, Ql) ~ ±g^{x, Ql) A6g \n{l/x) ± g^ix, Ql) ln(l/x), (66) 

where gi and g2 represent the two independent terms in the gluon parameterisation. The inter- 
play between the two terms allows for a large uncertainty at a; < 10~^ where the data constraint, 
from the Q^ dependence of F2{x,Q'^) at HERA, diminishes rapidly. 



A comparison of the fractional uncertainty for the present MSTW, CTEQ6.6 [87], Alekhin [139 
and NNPDFl.O ^1371] NLO gluon distributions is shown in Fig. [161 The fractional uncertainty for 
MSTW blows up very quickly at small x, whereas that for the Alekhin fit g ets bigger relatively 



slowly (the same would be qualitatively true of the HI |109l | and ZEUS jl40j fits), while that for 
CTEQ saturates, or can even decrease slightly, completely contrary to what is expected from 
the degree of constraint from data. These results are related to the input forms and values of 
Ql, and can be understood fairly easily. MRST/MSTW parameterise the starting distributions 
at Ql = 1 GeV^, and since we allow both positive and negative small-x contributions the un- 
certainty can be very large. We believe that this represents the true uncertainty at low x. The 
Alekhin, HI and ZEUS gluon distributions are input at a higher scale. They behave like x^'^ at 
small X, so the uncertainty is due to the uncertainty in one parameter, as shown in fl63|) . Recent 



CTEQ gluon distributions are input at Ql = 1.69 GeV^, and behave like x^ at small x where A is 



!0 

|18| 



large and positive, i.e. the input gluon is valence-like o This requires fine tuning: when evolving 



^^Despite the NNPDFl.O input parameterisation formally behaving as xg ^ x '^■'^ NNg{x), the neural netw ork 
parameterisation NNg(a;) should be sufficiently flexible that there is no dependence on the f actor x~° '^ |137 |. 



study of the input gluon parameterisation by CTEQ was made in Appendix D of Ref. [l26'|, but it did not 
address the impact of the input parameterisation on the uncertainties. 
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Figure 1 6: A comparison of t he fr actional uncertainty for the present MSTW, CTEQ6.6 87 



Alekhin [139| and NNPDFl.O [137|] NLO gluon distributions at = 5 GeV^. All uncertainty 



bands represent a 90% CL. limit. 



backwards from a steep gluon at higher scales a valence-like gluon only exists for a very narrow 
range of (if at all). In this case the small-x input gluon distribution is tiny, so there is a 
very small absolute error, as seen from flMl) . At higher all the uncertainty is due to evolution 
and driven by the higher x and more well-determined gluon distribution. Hence the very small 
X gluon is no more unce rtain than at x = 0. 001-0. Olo This does not seem to be realistic to 



us. The NNPDFl.O [l37l | uncertainty shown in Fig. [161 which is obtained using a very different 
approach for the input parameterisation and error propagation, will be discussed in Section [TH 

The above discussion of the input gluon parameterisation applies to all NLO and NNLO 
MRST/MSTW sets since 2001. However, in the MSTW 2008 analysis we have one more free 
parameter in the uncertainty determination, i.e. as well as 5^, 5gi and rjg we also have rjgi. The 
latter contributes to the uncertainty at intermediate x since its central value is much greater 
than that of rjg. The comparison of the uncertainty in the gluon distribution to that from the 
MRST 2001 and 2006 analyses is shown in Fig. [T71 As well as increased uncertainty at small x 
one sees that the "neck" at a; ~ 0.13 has disappeared in the new analysis. 

Note that the fractional uncertainty for the gluon distribution from the MRST 2006 NNLO 
set, shown in Fig. [TTT b). has an additional narrowing or "neck" at x = 0.7. This has a sim- 
ple explanation. The MRST 2004 [18] and 2006 [21] PDF^ts defined the input MS gluon 
parameterisation via a transformation from the DIS to the MS scheme, where the input gluon 



^^The same features exist for the "dynamical parton distributions" of Ref. |141i] . However, in that case the 
starting scale is Qq = 0.5 GeV^, and the PDFs do not fit data as well as the "standard" PDFs which are obtained 
from evolution at a higher starting scale. 
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(a) 



Gluon distribution at = 100 GeV^ 




(b) 



> 0.2 



Gluon distribution at = 100 GeV^ 




Figure 17: The fractional uncertainty on the gluon distribution at = 100 GeV^ (a) from the 



MSTW 2008 NLO fit compared to that from the MRST 2001 NLO fit M and (b) from the 



MSTW 2008 NNLO fit compared to that from the MRST 2006 NNLO fit [21 



distribution in the DIS scheme was given by the standard parameterisation (fTOl) . giving an en- 
hancement of the MS gluon at high x originating from the high-x quarks. Hence, in practice, 
at high X, i.e. above x ~ 0.4-0.5, the MRST 2006 gluon distribution was dominated by the 
quark distributions via the DIS MS transformation. This meant that for x > 0.6 the relative 
uncertainty on the gluon effectively became that on the quarks, which is much smaller, but does 
start diverging again at even higher clearly be seen in Fig. [T7l(b). As we will see in 

Section [121 the high- a; enhancement in the gluon distribution is disfavoured by Tevatron Run II 
jet data, and so the present MSTW 2008 fits revert to the standard parameterisation ffTOl) for 
the input gluon distribution directly in the MS scheme. Hence the uncertainty simply expands 
for the other PDFs. 

The LO gluon distribution has no tendency to turn over at small x, so the second term in 
(fTOjl is redundant. In this case we use the standard "MRS-type" parameterisation (160|) . but with 
all 4 parameters {6g, rjg, eg, 7p) free in determining the eigenvectors. 

Sea asymmetry, d — u 

The sea asymmetry d — u is parameterised at input as 

x{d- u){x, Ql) = A^a-u) x"'^"-^^ (1 - xf^^\l + 7(j_^) x + x'), (67) 

and the parameters contributing to the eigenvectors are the same as in the MRST 2001 fit. The 
distribution is constrained by data for 0.01 < x < 0.2. At small x the fact that we choose a 
simple power means that once the data for 0.01 < x < 0.1 cause the distribution to head to 
zero there is no scope for it to do otherwise at smaller x. We see no reason to suppose that, 
as a nonsinglet quantity, it should not go to zero as x goes to zero, due to considerations from 
Regge theory. At high x we choose the (1 — x) power to be rjs + 2, since d{x, Q^) — u{x, Q^) 
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Figure 18: The fractional uncertainty on the d — u distribution at = 100 GeV^ from the 



MSTW 2008 NLO fit compared to that from the MRST 2001 NLO fit [16 



is becoming very small, and we attempt to constrain u,d > 0. Applying this constraint on the 
power of (1 — x) we then need the term instead of the more usual y/x term to give sufficient 
flexibility at large x. The 3 eigenvector parameters control the normalisation, high x and low 
X. In fact, in some PDF sets one does then get negative central values of the d{x, Q"^) at the 
highest X, but at a level very small compared to the uncertainties. The uncertainties for d — u 
in the MSTW 2008 PDF set are slightly smaller than for MRST 2001, as seen in Fig. [TBI but 
this is due to the inclusion of new data and also smaller values of the tolerance. 

Light sea quarks 

The sea distribution S = 2u + 2d + s + s is parameterised by the standard "MRS-type" form, 

xS{x,Ql)= Asx^''{l-x)''^{l + ts^+^sx), (68) 

In this case there is no sum rule constraint, so in principle there are 5 free parameters. Again, 
only 3 can contribute to the eigenvectors due to very large correlations. We would expect these 
to be ^5, r]s and e^, but in practice 5s has too large correlations with the gluon parameters, so 
we take As instead. Hence, in this case, at small x the uncertainty at input is due to As-, and the 
uncertainty AS'(x, Qq) is simply proportional to S{x, Ql). However, at higher Q"^ the uncertainty 
is controlled by the gluon evolution which leads very quickly to a much larger uncertainty at 
very small x. This feature is shown in Fig. [TOT a). Perhaps we underestimate the sea uncertainty 
at small and very small x with this choice of free parameters. However, the uncertainty does 
follow the data constraint, i.e. ^2(0;, Q^) constrains the quark distributions down to a; = 10~^ 
quite precisely but only for a narrow range of small Q^. In Fig. [TW b) we compare the sea 
distribution at Q"^ = 100 GeV^ to that from the MRST 2001 analysis. The uncertainties are 
larger in the new analysis mainly due to the extra freedom in the strange quark parameterisation, 
which feeds into the uncertainties for u and d, and also due to allowing data set normalisations 
to go free; see Fig. [1] 
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Figure 19: The fractional uncertainty on the hght sea quark distribution S = 2u + 2d + s + s at 
Q2 = XOO GeV^ from the MSTW 2008 NLO fit compared to (a) the input sea quark distribution 
at = 1 GeV^ and (b) the sea quark distribution from the MRST 2001 NLO fit [lej. 

Strange quarks 

As mentioned before we now parameterise the strange quarks separately rather than assume 
that 

s{x, Ql) + six, Ql) = K [uix, Ql) + d{x, QD] , (69) 

with K ^ 0.4-0.5. We could use a completely free parameterisation, but the data only exist 
for X > 0.01 so this would lead to an enormous uncertainty for x < 0.01. This is essentially 



what is done for the CTEQ6.6 distributions [87|, but it is questionable whether this is realistic 



As a rough constraint, for example, Regge trajectory considerations suggest that all fiavour 
distributions have the same power as x — 0, and we certainly think that it is sensible to impose 
this. However, we apply more detailed constraints from perturbative QCD reasoning. The 
strange quark has some non-insignificant mass, and we assume that it is this which qualitatively 
leads to the suppression compared to up and down sea quarks. This assumption gives us a 
handle on the type of parameterisation to use. 

When charm and bottom quarks begin their evolution they evolve like massless quarks, but 
always lag behind. This leads to some suppression at all x for finite Q^. But at small x this 
suppression is roughly just a normalisation suppression, i.e. evolution makes the small-x powers 
Si the same. Indeed, other than the overall normalisation suppression, further effects are mainly 
at high X. This reasoning leads to us choosing only the normalisation and the high-x power 
ry+ as free parameters for s{x,Ql) + s{x,Ql). Indeed, effectively we choose a parameterisation 
of the form 

xs{x,Ql) + xs{x,Ql) = A+ x^''- {1 - xy+ {1 + es + Is x) 

= A{l-x)^xS{x,Ql), (70) 

where A = A^/As, fj = t]^ — rjs and the small-x power is fixed to 6s- Introducing distinct e and 
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MSTW 2008 NNLO (90% C.L.) 



^ at = 1 GeV^ 
% ^atQ^ = 15GeV^ 




(b) 



(s+s)/(u+d) distribution at = 5 GeV^ 




CTEQ6.6NLO 
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Figure 20: (a) s + s at = 1 GeV^ compared to c + c at = 15 GeV^, both as a fraction of 
u + d. (b) {s + s)/{u + d) at = 5 GeV^ from the current NLO analysis compared to the same 
quantity from the MRST 2001 [ij and CTEQ6.6 analyses. 



7 parameters do not improve the fit quality and lead to instability in eigenvectors, so these are 
also fixed to the same values as the total sea. Therefore, apart from the normalisation Aj^ and 
the high- a; power all other parameters in the s + s parameterisation are fixed to be the same 
as the total sea S. 

The results of the fit will be discussed in much more detail in Section [HI but at input there is 
about a 35% normalisation suppression of the strange distribution compared to the average of 
u and d at Ql and some additional high-x suppression. Hence the suppression at the borderline 
nonperturbative scale Ql = 1 GeV^ is now ~ 0.3 which is the value used in hadronisation 
models (i.e. the probability to generate ss pairs compared to uu or dd pairs). This is shown 
in Fig. [20](a) along with the suppression of c + c at = 15 GeV^, i.e. evolved through ~ 7-8 
times the mass scale. This charm suppression is very similar to that for s + s at = 1 GeV^, 
implying a suppression due to evolution from slightly more than 0.1 GeV^, similar to m^. We 
would not of course expect exact correspondence, especially since the strange evolution is in a 
nonperturbative regime; however, the comparison is surprisingly good except that c + c is more 
suppressed at x ~ 0.1 (which is the implication for s + s from recent HERMES data on 
production 142| ). This gives us confidence that the difference between the total sea and the 



strange distribution is mainly due to mass suppression, and that our limited parameterisation is 
well motivated. In Fig. [2U](b) we compare the ratio {s + s)/{u + d) at = 5 GeV^ obtained from 
the present analysis to the same quantity from the MRST 2001 16|] and CTEQ6.6 (87| analyses. 
The MRST 2001 ratio is slightly larger and has much smaller uncertainties, since the input was 
taken in the fixed form (!69|) with k = 0.5. The CTEQ6.6 ratio in Fig. [20] is much larger for 
X < 0.01, and has larger uncertainties, due to the more flexible parameterisation in a region 
where there are no constraining data, but is surprisingly also larger in the region 0.01 ^ x < 0.2 
where the dimuon data should provide a constraint. 
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The strange asymmetry is not at all well-constrained and the parameterisation, 

xsy{x, QI) = xs{x, QI) — xs{x, QI) = x^^ (1 — x)^' (1 — x/xq), (71) 

is the simplest we can think of with the correct features. The number sum rule of zero strangeness 
determines Xq, and hence it is a function of the other parameters. In practice the data are 
not sufficient to constrain and 5_ independently, so we fix the latter to 6- = 0.2. Eigenvector 
20 (at NLO and NNLO) consists mainly of hence its peculiar nature shown in Fig. [6l There 
is little constraint from data, but the allowed values of the parameter are limited. 



Summary 

In summary, we parameterise with simple forms which do not allow bumps or shoulders in 
general. Some distributions are restricted by theory assumptions, i.e. the distributions stay 
extremely small {d — u) or behave like mass-suppressed quarks (strange). We beheve that 
we obtain a good representation of the possible uncertainty in regions where there is a data 
constraint, since there is an enormous amount of constraining data. We concede that we may 
underestimate the uncertainty in regions of extrapolation, e.g. the small-x valence quarks, d — u 
and strange distributions. However, for some of these the uncertainty is essentially unlimited 
without making some theory assumptions. 



6.6 Fit to a reduced dataset 



As part of the first proceedings of the HERA-LHC workshop [71|, a "benchmark" PDF fit was 
proposed by S. Alekhin consisting of a limited number of DIS data sets with fairly conservative 
cuts of > 9 GeV^ and W'^ > 15 GeV^. The original aim was to provide a test of different fitting 
codes starting from a common input and theoretical assumptions and fitting to the same data. 



analogous to the benchmark tables for PDF evolution provided by G. Salam and A. Vogt |70l. |71 



The benchmark fit was also carried out by R. Thorne [71| using the MRST fitting code, albeit 
with some differences in the data sets used and the treatment of correlated systematic errors, but 
good agreement was found with the benchmark PDFs of S. Alekhin. The benchmark PDFs were 
also compared to those from the MRST 2001 NLO global fit, and significant discrepancies were 
found, well outside the error bands, suggesting some inconsistency of the data sets included in 
the global fit compared to those in the benchmark fit, and an inadequacy of the method for error 



propagation [TjJ. However, the benchmark PDF fit of Alekhin/Thorne used a more restrictive 
parameterisation for the input gluon and sea quark distributions than that used in the global 
fits. Also, the one-sigma PDF uncertainties for the benchmark fit were defined by Ax^ = 1, 
which is clearly inadequate. The benchmark fit was repeated using the M STW 2008 analysis 
framework as part of a contribution to the second HERA-LHC proceedings jl43| . Although the 
discrepancy between the benchmark and global PDFs was reduced compared to the previous 



study [TjJ, there were still sizeable discrepancies in the gluon at low-x and in the down valence 
distribution. It was realised that some of these discrepancies were due to different values of as 
used and due to the unnecessary assumption that d = u in the benchmark fit. It is therefore 
interesting to repeat the benchmark fit here, but fixing 0:5 to the value in the MSTW 2008 
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Data set 


Reduced fit 


Lriouai nt 




i/O / iO( 


1 77 / T c;7 


NMC /ip F2 [33] 


72 / 67 


75 / 67 


NMC [103] 


59 / 73 


65 / 73 


HI 96-97 e+j9 NC [109] 


46 / 80 


47 / 86 


ZEUS 96-97 e+p [112] 


155 / 206 


163 / 206 


All data sets 


458 / 589 


526 / 589 



Table 6: x^/^pts. for the data sets included in the reduced fit compared to the corresponding 
values in the global fit. 



global fit and allowing for d ^ u. In this way, we can investigate potential data set inconsistency 
and the (in) adequacy of the error propagation by fitting to a reduced number of data sets, but 
without any bias due to unnecessarily different input assumptions compared to the global fit. 
An a nalo gous investigation has recently been carried out by the NNPDF group in Section 5.5 of 



Ref. 137 



We use the same data sets as specified in the MSTW version of the HERA-LHC bench mark 
fit that is, we fit BCDMS data on (sij, NMC data on Q, NMC data on Fi / |l03 



HI 96-97 e+p NC data on a (loj and ZEUS 96-97 e+p NC data on d (nj, with cuts of > 9 



GeV^ and > 15 GeV^. The heavy fiavour treatment is the same as in the global fit, that is, 
using the GM-VFNS with rric = 1.4 GeV and rrib = 4.75 GeV. The strong coupling is fixed at 
the value obtained in the NLO global fit. The input PDF parameterisation is taken to be the 
same as in the global fit, with the following exceptions. None of the data sets included constrain 
strangeness, and indeed, negative strange quark distributions are obtained if allowed to go free. 
Therefore, we fix the input Sy = s — s and {s + s)/{u + d) to be the same as in the global fit. 

In Table |6] we compare the values for the reduced fit to the values for the same data 
points in the global fit. The deterioration in the quality of the fit (total A^^ = 68 ^ 1) with 
the inclusion of more data sets is evident, especially for the BCDMS F2 data. Both the reduced 
fit and the global fit give an acceptable description of all data sets included in the reduced fit, 
and so any error propagation based on Ax^ = 1 is clearly inadequate. 

Application of the dynamic tolerance procedure described above to the reduced fit gives 
values shown in Fig. [211 Typical values for the tolerance are T ~ 3 for 68% C.L. uncertainties 
and T ~ 5 for 90% C.L. uncertainties, that is, slightly smaller than the typical values in the 
global fit, and with less variation between different eigenvector directions. We note that the 
departure from ideal quadratic behaviour (T = t) for higher eigenvector numbers is more severe 
for the reduced fit than the global fit, since the PDF parameters are much less constrained in 
the reduced fit and hence there are more relatively "fiat" directions in eigenvector space. 

In Fig. [22] we compare the and d^ distributions at = 20 GeV^ from the reduced and 
global fits. The dashed lines indicate the uncertainty bands obtained with T = 1, while the outer 
uncertainty bands are obtained using the dynamic tolerance shown in Fig. [HI For the valence 
quarks one can see that there are differences in detail, particularly for d^. Using uncertainties 
determined by Ax^ = 1 there are still significant discrepancies in places, e.g. d^j aX x = 0.05, 



64 



20 
15 
10 
5 

-5 
-10 
-15 
-20 



MSTW 2008 NLO (reduced) PDF fit 



Q. a 

■ « o ■ 

1- CD 



D 

■ G- 



a 



.JO) O. . 

T- o 
X m 



■5^ <^ 



a 
o 
m 



UJ 



a 
a 
ffi 



o 
o 



..Or.. 

u> 
s 

Q 
O 

m 



■1^ 

-Q. - - Q- 

0) o 
T- m 



+\100 (CTEQ) 
+\|50 (MRST) 



...j -VSO (MRST) 
»(lOO (CTEQ) 



7 8 9 10 11 12 13 14 15 16 

Eigenvector number 



Figure 21: Tolerance for each eigenvector in the reduced fit. 

but using the dynamic tolerance the two fits are consistent, though the global fit leads to much 
smaller uncertainties, as expected. In Fig. [231(a) we see the curves for the u distribution, where 
again there is good compatibility, though the lack of Drell-Yan and neutrino DIS data in the 
benchmark fit leads to a huge uncertainty at high x. At a; > 0.5 there is some discrepancy. This 
is related to the situation for the gluon, seen in Fig. [231(b). At small x the two gluons agree 
well, again with that from the reduced fit having larger uncertainty. At high x there is some 
disagreement. Using as{M'^) ~ 0.120 the extremely good fit to high-x BCDMS data, seen in 
Table [6], requires a large high- a: gluon in order to obtain the sufficiently flat shape in preferred 
by the data. This produces a pull in clear contradiction to other data, and which is overwhelmed 
in a global fit, but survives in the reduced fit. 

Hence, we can see that the "benchmark" type of fit to a reduced number of data sets proposed 
in Ref. 71|is a useful means of comparing different fitting approaches, as demonstrated in 
Refs. 7l[ . [l43[ |. However, the central values of the PDFs obtained and the uncertainties obtained 
using the Hessian method with the textbook Ax^ = 1 should not be taken seriously — hence our 
decision not to make these sets publicly available. As shown above, comparison of the benchmark 
PDFs with those obtained from more complete fits can be useful in suggesting the reliability 
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(a) (b) 



Up valence distribution at = 20 GeV^ Down valence distribution at = 20 GeV^ 




Figure 22: Comparison of the (a) up and (b) down valence quark distributions from the global 
and reduced NLO fits. 
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(a) 



(b) 



Up antiquark distribution at Q = 20 GeV^ 

MSTW 2008 NLO (68% C.L.) 
Fit to reduced dataset 




lUISTW 2008 NLO (68% C.L.) 
Fit to reduced dataset 




Up antiquarl< distribution at = 20 GeV^ 



Gluon distribution at = 20 GeV^ 



IVISTW 2008 NLO (68% C.L.) 
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IWSTW 2008 NLO (68% C.L.) 
Fit to reduced dataset 




Figure 23: Comparison of (a) the u distribution and (b) the gluon distribution from the global 
and reduced NLO fits. 
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of more elaborate methods of determining the PDF uncertainty, giving good support for the 
"dynamic tolerance" approach introduced in this paper. Benchmark fits can also be useful for 
highhghting data sets which are inconsistent with the rest of the data in a global fit. Indeed, 
it is shown that there is clear tension between the BCDMS data and the rest of the data, 
which is a concern, although the degree of compatibility becomes better when the cut of > 9 
GeV^ is lowered to our standard choice of > 2 GeV^. 



68 



7 Description of inclusive structure function data 



7.1 HERA structure functions 

We show the fit to a selection of the HERA neutral-current data for the reduced cross section, 

Q') = F^ix, Q') - , ^ ,f Fl{x, Q'), (72) 
1 + (1 - y^) 

in Fig. [231 This illustrates the general features, and the quality of the fit can be judged from 
the values in Table [2] in Section [5l At each order the quality of the fit is very similar to 
previous MRST fits. A number of general features may be observed by examination of Fig. [241 
It is clear that the LO fit is not able to generate sufficient evolution to give a genuinely good 
fit, being fiatter at both low and high in each x bin. This is due to the presence of large 
positive contributions to the splitting function Pgg beyond LO, and even the very large value of 
«5 in the LO fit cannot fully compensate. We also note that the turnover of 5"(x, Q^) at high 
Q'^ (equivalent to high y) in each x bin is largest at LO, particularly at small x. As we will see 
later in Section [T3] this is because Fl{x,Q'^) is very large at LO. The NLO and NNLO fits are 
very similar to each other. The main difference is that the NNLO curves are rather steeper at 
low X and Q^, being driven by an extra ln(l/a;) factor in the NNLO Pgg. Although it is difficult 
to see fror n th e plot this does actually lead to a worse fit at NNLO for the HI 97 minimum bias 



e+p data |l09l |. see Table [21 the data preferring a slightly fiatter shape. Overall, the quality of 
the fit is very good. The only fiaw is a slight trend for the evolution to not persist quite strongly 
enough at high Q'^ for x ~ 0.001-0.01, although there is quite a large spread between the HI 
and ZEUS data at high in these x bins. This might possibly be a sign of further required 
corrections to the theory of some sort, and indeed an improvement due to small-x resummation 



was shown to be possible in Ref. |144 



The only HERA structure function data completely new to our analysis are the HI MB 99 
e'^p data [Sli — a handful of points which probe the highest y values and hence display sensitivity 
to Fl{x, Q'^). In Fig. [25lwe show our results in bins of to best compare to these specific data. 
The dashed lines show the one-sigma PDF uncertainty bands. We noted in Ref. 28[ that while 



the turnover at small x (high y) did not occur with our central NLO fit, it did occur with the 
central NNLO fit due to the increase of Fi{x,Q'^) from the NNLO coefficient function. It is 
striking that the central value of the NNLO result no longer fully turns over at low and small 
x. However, from Table [21 one can see that the fit to the HI MB 99 e^p data is still a little better 
at NNLO, but also that it is not very bad at LO or NLO. The data, while being indicative, are 
not in practice very constraining. It is also clear from the plot that the full turnover is allowed 
within the uncertainties of the PDFs. We will discuss more direct results about, and relating 
to, Fl{x,Q'^) in Section [131 



Charged-current data from e~^p scattering 35|, |36|] are also included in the fit. They are fit 



perfectly well, as seen in Table [21 but at present do not provide a strong constraint on the fit. 
Much more precise data are expected from HERA II analyses. 
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Figure 24: The reduced cross section a{x, Q"^) for a selection of HERA data and the structure 
function ^2(0;, Q^) for a selection of fixed-target data. The curves represent a(a;, Q^), evaluated 
at the HERA centre-of-mass energy, and so for the high x and low values corresponding to 
fixed-target data, y is negligible and the curves effectively represent the required -^2(0;, Q^). The 
normalisations of data sets are those in the best fit at NNLO, found in Table [31 These are very 
similar to those at NLO, but at LO the HERA data would be 0.5-1% lower in general. 
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7.2 Fixed-target neutral-current structure functions 

The fit to a selection of fixed-target data is also shown in Fig. |24l where in this case the data 
and the curve (to an extremely good approximation) represent the structure function F2{x, Q^). 
These data have appeared in numerous previous MRST analyses, and at NLO and NNLO the 
fit is of the same quality as previously. These data are one of the main constraints on the 
value of as'(Ml), though it is correlated to the high-x gluon. Indeed, the value of as{M^) 
at NNLO is slightly smaller than in our most recent analyses f2l|], and so the high-x gluon is 
correspondingly smaller. Perhaps the most striking new feature is that at LO the value of 05 
has increased sufficiently that the fit quality is much closer to that at NLO and NNLO, the 
size of the coupling largely making up for the increased evolution at higher orders in the high-x 
regime. Indeed, at x = 0.35 we see that at very low (below our W'^ cut on data) the LO 
curve in Fig. [2l] is falling most quickly — a feature not seen previously. 



7.3 Neutrino structure functions 



There have been significant changes in both the data to which we fit, and the details of the 



procedure used, for neutrino structure functions. The NuTeV data [37 



on 



and X-F3 = (xFg + xF^ )/2, where N = {n + p)/2, replace the previous CCFR data |145 



These two experiments both use an iron target, corrected to an isoscalar target, and cover much 
the same kinematic range, but the NuTeV data are more precise. The two sets agree fairly 
well, except in the very high-x (i.e. x > 0.5) region, where the main source of the discrepancy 
is the calibration of the magnetic field map of the muon spectrometer. In practice we find the 
high-x NuTeV data very difficult to fit. At high-x the predictions are mainly determined by the 
valence up quark distribution, which is very well constrained by the neutral-current structure 
function data. Given the degree of experimental uncertainty in this region, we do not fit to 
neutrino data with x > 0.5. The main information from these data comes from x < 0.3, where 
we are sensitive to the details of the different valence and sea quark contributions, which are 
weighted differently for neutrino structure function data compared to neutral-current data. We 



also include the recent CHORUS F2 and xF^ data |38|] which cover a similar range in x but at 
slightly lower on average, and are taken using a lead target. These are completely compatible 
with the NuTeV data for x < 0.5, but seem to lie a little lower at the highest common x value of 
0.65. We apply the same cut to the CHORUS data of x > 0.5 for consistency. As mentioned in 
Section [5] we apply a cut of W"^ > 25 GeV^ on the xF^ data to guard against potentially larger 
(than F2) higher- twist contributions. 

In the analysis of these data, we have adopted a more sophisticated treatment of nu- 
clear corrections compared to previous fits. We apply the nuclear corrections Rf, defined as 
f^{x,Q'^) = Rf{x,Q'^, A)f{x,Q'^), separat ely fo r each parton flavo ur / using the results of a 



NLO fit by de Florian and Sassot (nDS) [mJ, using the GRV98 [l48j free proton PDFs as 
input; see Fig. [261 Here, f"^ are defined to be the parton distributions of a proton bound in 
a nucleus of mass nu mbe r A. In the figure we compare to a NLO analysis by Hirai, Kumano 
and Nagai (HKN07) using the MRST98 freeproton PDFs as input, and to the simple 



"RATFE" fiavour- and Q^-independent corrections [12], applied in previous MRST fits. Since 
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Figure 26: Nuclear corrections for iron {A = 56, Z = 26), at = 20 GeV^. The solid lines, and 
the associ ated one-sigma uncertainty bands, are from an analysis by Hirai, Kumano and Nagai 
(HKN07) |l46 |. The dashed lines are from an analysis by de Florian and Sassot (nDS) jl47 |. and 
the dotted lines are multiplied by the correction factor, fl75]) . determined from our NLO global 
fit. We take the nDS nuclear corrections for charm (and bottom) quarks to be the same as those 
for strange quarks. The flavour- and Q^-independent "RATFE(x)' 
the dot-dashed lines, were applied in previous MRST flts. 
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the nDS analysis [l47l | neglects heavy flavours, we take the nuclear corrections for charm (and 
bottom) quarks to be the same as those for strange quarks. The nDS and HKN07 nuclear correc- 
tions are similar for the up and down quarks, which are most important for the neutrino structure 
functions, but differ significantly for the less important gluon and other quark flavours An 
alternative approach consists of extracting nuclear PDFs from data on a single nuclear target 
(iron) [l5n| . 



The nDS nuclear corrections do not have an associated uncertainty. We assume that the 
uncertainty is likely to be of order a few percent, and account for it by multiplying the nuclear 
corrections Rf{x, Q^, A) by a flavour-independent modification function of the form: 

NUCMOD(x) + + ^ (73) 

|(1 + 0.03 ri)[l + 0.015 rs ln2(x/x^)] : x > Xm 

where Xm = exp(— 2.5) is chosen to be roughly in the middle (logarithmically) of the x range 
spanned by the neutrino data. The parameters ri, r2 and are allowed to go free when we 
obtain our fits. This function allows, for example, a little more variation than the uncertainties 
in the correction factors quoted in HKN07 [l40|. The extent to which the central values of the 
nDS nuclear corrections are modified by fl73l) is shown in Fig. [26j 2-3% is typical, and we regard 
this as perfectly acceptable. 

The quality of the NNLO fit to F2 data is shown in Fig. [571 where the data which fail to 
satisfy the W"^ cut and the x < 0.5 cut are also shown. The fit quahty is perfectly good, as 
seen in Table [2], and the NuTeV and CHORUS data agree very well in the region of overlap, the 
latter lying a little lower in Q^. For x > 0.5 there is a slight tendency for CHORUS data to lie 
above the prediction, which is dominated by the quark determined by neutral-current DIS 
data. The NuTeV data, however, lie consistently well above the prediction. 

The quality of the NNLO and LO fits to xF^ data is shown in Figs. [28] and [29], respectively. 
One obvious result of the LO fit is that the slope of the theory is consistently too fiat, except 
for X > 0.275, and the fit is correspondingly poor. This failure is largely cured in the NNLO fit, 
though perhaps not completely at the lowest x values. Indeed, it is clear from the values in 
Table [2] that the fit to 0:^3 data, for both experiments, improves significantly from LO to NLO, 
and then even more at NNLO. Since 

xF, ^ i (xF3^^ + = C3 ® 5^(g. - g.), (74) 

i 

this is a particularly clean test of perturbative QCD (and as) with no complications from mixing 
in parton evolution — once the single parton combination is determined at one scale we have 
a complete prediction (in terms of as) at higher Q^. In this case the improvements due to 
the change in shape with of the structure function at increasing order are clearly seen. As 
seen from Table [21 the values for the 45 NuTeV xF^ points improve from 62 — 40 — ^ 34 
in going from LO— i>NLO— s>NNLO, while the values for the 33 CHORUS xF^ points improve 
from 44 — * 31 — > 26. It is not certain that the improvement has converged at NNLO, and 



20The EPS09 [m^ corrections are in good agreement with those from nDS [l47l | for valence and sea quarks 
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Figure 27: 
CHORUS 
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The quality of the NNLO fit to F2 = F^'^ = F^'^ data from NuTeV [37] and 
The data which fail to satisfy the W"^ cut and the a; < 0.5 cut are also shown. 
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Figure 28: The quality of the NNLO fit to xF^ = {xF^^ + xF^^)/2 data from NuTeV [37[ and 
CHORUS Ssj. The data which fail to satisfy the W'^ cut and the x < 0.5 cut are also shown. 
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Ratio of NuTeV and CHORUS xF, data to MSTW 2008 LO PDF fit 
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Figure 29: 
CHORUS 



= {xF^^ + xF^^) /2 data from NuTeV |37| and 



The quality of the LO fit to 0:^3 

381]. The data which fail to satisfy the W'^ cut and the x < 0.5 cut are also shown 
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there is perhaps a httle scope for higher-order corrections to improve matters further. This is 
quite possibly the case for many quantities, but the mixing between parton flavours obscures the 
interpretation for most. It is also interesting that these data favour fairly high values of a5(M|), 
contra ry to the claim often made about the results obtained from fitting structure function data 



alone [1391 . Il5ll . Il52l |. The trend of the xF^ data at x > 0.5 to lie above the theory predictions 
is the same as for F2. In Figs. [28] and [29] we also show the data below our higher (than for F2) 
W"^ > 25 GeV^ cut. We see that for x < 0.35 these data do indeed lie systematically below the 
NNLO result, implying a large negative higher-twist contribution, which may also very slightly 
influence data above our cut. We leave a systematic analysis of such higher- twist contributions 
to a future study. 
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8 Neutrino dimuon production: constraints on strangeness 



Our new analysis allows much more flexibility for the strange quark and antistrange quark 
distributions. This is both for the total and the two separate contributions. Protons have no 
valence strange quarks, that is, 



dx [s(x,g2) - s(x,g2)] = 0. 



(75) 



But this does not necessarily mean that s{x, Q"^) = s{x, Q"^) for all x. For example, in a meson 
cloud model, the proton can fluctuate into a A{uds) and a or K*{us), where the s in the 
K or K* carries a larger momentum fraction of its parent hadron than the s quark in the A. 
Nevertheless, the assumption in earlier PDF fits has been that the input distributions satisfy 



(u + d), 



(76) 



with a constant k ~ 0.4-0.5, justified by a simplified analysis of dimuon data by the CCFR 
experiment [4Ci]. Updated CCFR and NuTeV dimuon cross sections are now available, therefore 
we include these data in the global fit to constrain s and s separately; see, for example, the 
process Ufj^N fi^fi'X proceed ing via t he subprocess W~^s —>■ c shown in Fig. [301 Previous 
studies have been made in Refs. 153- 156 . 




Figure 30: Dimuon production in neutrino-nucleon DIS from charm production via scattering 
of the boson off a strange or down quark from a nucleon in the heavy nuclear target N. 
Diagram taken from Ref. 39 . 
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8.1 Description of dimuon cross sections 

can be related to the underlying 



The CCFR and NuTeV dimuon cross sections, v^N 



process, v^N 



^ cX, using 

dV ^ 
dxdy ^ 



dV 
dx dy 



(77) 



and similarly fo r th e antineutrino data. Here, Be = 0.099 ± 0.012 is the charm semileptonic 
branching ratio 1531]. The kinematic acceptance corre ction A accounts for the 5 GeV cut on the 
energy of the ch arm d ecay muon and was calculated 153l | for each bin usin g th e DISCO Monte 
Carlo simulation 157 ] assuming a charm mass of 1.4 GeV and Collins-Spiller |l58 | fragmentation 
parameter e = 0.6. Allowing for the CKM mixing, the cross section has the 0{a^) form: 



A 



A 



(7J 



where the rescaling va riabl e ^ = x{l + ml/Q'^), and where Z = 23.403 and A = 49.618 for the 



CCFR/NuTeV target |153l |. Thus, the and dimuon cross sections are direct measures of 



s and s, respectively. We have assumed isospin symmetry, dT' = = u, and, as usual, have 
taken the renormalisation and factorisation scales to be fiR = = Q- Finally, just as fo r the 
inclusive neutrino data, we apply the nuclear corrections of de Florian and Sassot (nDS) jl47 
for each parton flavour, and hence relate the nuclear PDFs, /^(x, Q^), to the proton PDFs, 



f{x, Q^); see Section rr3l As evident from the difference between the nDS |l47l | and HKN07 |l46 
parameterisations for nuclear corrections to strange quarks shown in Fig. [26l the uncertainty is 



likely to be a few percent. However, this uncertainty is smaller than the experimental uncertainty 
on s + s at the relatively low scales where the data are probed. Additionally, there is no evidence 
for any systematic problem with shape in the flt quality. 

It is convenient to parameterise s ± s, rather than s and s separately, as explained in Section 
[6l Recall that at the input scale Qq = 1 GeV^ we take 



x{s + s) = A+ x^^l-xf+il + esV^ + ^sx), 
x{s -s) = A_ x^- (1 - xf'{\ - x/xq), 



(79) 
(80) 



where the parameters 5si and 75 are in common with the parameterisation of the total light 
quark sea, S, see ([8]). The flnal factor in flHOl) is to ensure that the number sum rule fl75|) is 
satisfled, and implies that (s — s) changes sign at a; = Xq. 

The description of the CCFR and NuTeV neutrino and antineutrino dimuon data given by the 
NNLO flt is shown in Figs. [3T], [321 [33] and [311 Clearly the quality of the flt is very good. There is 
perhaps a slight tendency for the CCFR antineutrino data to lie a little above the curves and the 
NuTeV antineutrino data to lie a little below the curves, suggesting that the more recent NuTeV 
data favour a lower s(x, Q"^) distribution in the region of the data and are mainly responsible 
for the asymmetry. The results are similar at LO, NLO and NNLO, as suggested by the quality 
of the flts shown in Table [21 Note that the branching ratio Bf. = 0.099 is flxed in the flts. If 
allowed to go free there is a distinct tendency for it to choose a high value. This is the same sort 
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Figure 31: The quality of the NNLO fit to the CCFR neutrino- initiated dimuon production. 
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;ure 32: The quality of the NNLO fit to the CCFR antineutrino-initiated dimuon production. 
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;ure 33: The quality of the NNLO fit to the NuTeV neutrino-initiated dimuon production. 
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;ure 34: The quality of the NNLO fit to the NuTeV antineutrino-initiated dimuon production. 
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Figure 35: Ratio of the strange to the non-strange distributions at the input scale Ql = 1 GeV^ 
and the associated one-sigma uncertainty band for both (a) NLO and (b) NNLO distributions. 



of effect as data set normahsations sometimes floating to one end of the allowed uncertainty. 
The combination 4:/9{u + u) + l/9{d + d + s + s) is very well constrained by structure function 
data. The fit to this can be obtained with the minimum momentum of partons if there is a 
large proportion of up and antiup quarks, so in the absence of higher order corrections the fit 
tends to choose this scenario if g iven the choice, and hence push the strange proportion down. 



Indeed it was recently noted |159l | that if the strange distribution is left completely free in a fit to 
structure function data alone it settles on a marginally negative value (we have also observed this 
phenomenon), albeit with very large uncertainty. We remove this tendency from the fit by fixing 
the branching ratio at = 0.099 due to complete correlation with the (s + s) normalisation. 

8.2 Impact on strange and antistrange distributions 

In Fig. [35] we show the ratio of the strange sea to the non-strange sea, 

s{x,Q^) + s{x,Q^) 



k{x, Q 



uix,Q^)+d{x,Q^) 



at the input scale Qq = 1 GeV^ at both NLO and NNLO. At small x this suppression factor 
is just a little below the previous default values of 0.5 at NLO and 0.44 at NNLO. However, 
there is a strong suggestion of additional suppression of the strange sea at large x relative to 
the non-strange sea. This is as one might expect if the suppression is due to the mass of the 
strange quark since at high x one is nearer to the production threshold. As seen in Section 16.51 
and illustrated in Fig. [20ra). the suppression of the strange distribution is consistent with the 
expectation from considering the strange quark to have a mass ~ 0.1 GeV^. In Fig. [361 we 
also show the difference between s and s at the input scale. Note that s and s are constrained 
by the dimuon data for 0.01 ^ x < 0.2, the region where the central value of the asymmetry is 
non-zero. However, note that the asymmetry is non-zero just outside the one-sigma uncertainty 
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Figure 36: The strange asymmetry, xs — xs, at the input scale Ql = 1 GeV^ and the associated 
one-sigma uncertainty band for both (a) NLO and (b) NNLO distributions. 
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Table 7: Integrated strangeness momentum asymmetry and integrated strangeness fraction for 
g2 = 1, 10, 100 GeV^. Note that J^dx x{s - s) ^ 0.007 (for typical ~ 20 GeV^) is required 
to bring the NuTeV sin^ 9w to the world average value. The one-sigma PDF uncertainties are 
given. 



band, and it is certainly consistent with zero within the 90% C.L. uncertainty band. As with the 
total strange distribution the plots for the NLO and NNLO analyses are very similar. Indeed, 
this is quantified in Table [7| which shows the integrated strangeness momentum asymmetry, 

[ dxx[s{x,Q^)-s{x,Q^)], (82) 
Jo 

and the relative strength of the strange quark sea 

J^dxx[s{x,Q^) + s{x,Q^)] ^g^^ 
/g dx X [u{x, Q"^) + J(x, Q^)] ' 

for Q2 = 1^ 10, 100 GeV^ for the LO, NLO and NNLO parton sets. Notice that in both 
cases there is good perturbative convergence, though the LO strangeness fraction is a httle low 
compared to the higher orders. This was not guaranteed to be the case, and relies on a careful 
definition of the GM-VFNS for charged-current structure functions, as outlined in Section HI 
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The uncertainty on the s and s distributions has increased significantly due to this new 
treatment, since previously it was simply the same (proportionally) as for the much more accu- 
rately determined u and d distributions. This increase is illustrated clearly in Fig. [201(b). and 
is one of the main reasons for the increase in uncertainty for the total sea seen in Fig. [191 The 
increased strange uncertainty also leads to a slight increase in the uncertainty for the u and 
d distributions. From constraints on the different weighting in the sum of sea quarks, mainly 
on the charge- weighted combination in neutral-current DIS, there is an anti-correlation between 
strange and u, d. However, this is a fairly small effect, and is difficult to disentangle from the 
consequences of other changes in our procedure (e.g. accounting for data set normalisation un- 
certainties). However, as well as an understanding of the uncertainties we now have much more 
confidence in our central values for the strange and antistrange distributions. The precise val- 
ues are important for current experiments since the strange quark and antiquark distributions 
make significant contributions to and W~ production at hadron colliders, particularly at the 
LHC, through cs and cs partonic fusion. The asymmetry is driven mainly by the difference 
between ud and ud fusion, but also has some sensitivity to the strange quark asymmetry. It is 
worth noting that the associated production processes pp W~ + c + X and pp + c + X 

are a promising way to measure the strange quark distributions, where the LO partonic subpro- 



cesses are gs ^ W c or qg — » W ^c [l56|- Indeed, the Wc production process has recently been 
measured at the Tevatron IGOl . 



8.3 Implications for the NuTeV sin^^^ anomaly 

The determination of the momentum asymmetry ^^dx x{s — s) of the s and s distributions 

reported by 



has imp orta nt implications 
NuTeV 



161 



162 



for the anomaly in the measurement of sin 9w 
The NuTeV extraction of sin^ Ow from 



R- 



- sin^ 9 



fi-X) -(r{u^N ^ fi+X) 2 ^ ^^^^ 

is about ~ three-sigma above the global average, and was at first thought to hint at new physics. 
However, there are two effects sensitive to parton distributions which should first be included: 
isospin violation [u^ ^ ci", (dF ^ u") and the strange sea asymmetry (s ^ s). Isospin violation 
gives 



R- 



+ (1 



^ sin^ Ow 



dx X [« 



<) 



2 /q da; x{u„ + d^) 

In fact, isospin violation is automatically generated by QED corrections to parton evolution, and 
was found by MRST 19|] to remove slightly less than half of the total discrepancy in the NuTeV 
value of sin^ 6w if evolution of the photon distribution (in the proton) from current quark mass 
scales was applied, while for evolution from constituent quark masses a reduction of about a 
quarter was obtained. The second effect is the possible existence of a momentum asymmetry 
j^dx x{s — s) of the strange sea. It gives 

7 . . /o^dx 



R- 



— sin 6w — (1 



sin 6y/) 



x[s 



/q dx x{uv + dv) 



86) 



87 



A positive val ue of J^dx x{s — s) ~ 0.007 1631 . Il64 | would bring the NuTeV sin^ 6w to the 



world average 132| . Indeed, we see that the analysis of the dimuon data gives a preference for 



a positive strange momentum asymmetry, which goes some way to further reduce the NuTeV 
sin^ 9w anomaly. A reduction in the sin^ 9w anomaly of about one-sigma is likely, but a reduction 
of up to two-sigma (or down to no effect) is within the range of the PDF uncertainties for the 
strange momentum asymmetry shown in Table [71 
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9 Description of Ff and data 
9.1 Comparison to HERA data 



The new Fn*^ data included in the MSTW fit 41-47 extend the data included in earlier MRST 



fits 4l|, |45j. In Fig. [37] we show the comparisons of our results at LO, NLO and NNLO to the 
charm structure function F2'^{x, Q"^) for all available HI 

liil-Q and ZEUS jii-Q data. As one 
can see from the figure, and from the values of iii Table [21 the quality of the fit is good. At LO 
the fit undershoots badly for = 2 GeV^, the absence of the 1/x term in the NLO coefficient 
function being most important in this region. However, the quick evolution at LO due to the 
large gluon and very large as quickly compensates and F2'^{x, Q^) evolves most quickly at LO. 
The NNLO fit quality is best. This is clearly due to the improvement compared to NLO at the 
lowest X for = 2 GeV^. At higher the NNLO F2'^{x,Q'^) evolves most slowly, a feature 



that was highlighted as a consequence of heavy flavour at NNLO in Ref. [72|. This is most 
easily understood in the GM-VFNS, where the positive corrections to the coefficient functions 
at low Q^, and the negative contributions to the heavy flavour distributions from matching 
conditions influencing high Q^, both have the effect of flattening the slope. However, it has 
recently also been noticed in an approximate FFNS at NNLO (inclusion of threshold corrections 
and factorisation/renormalisation scale dependent terms) j95|. This slower evolution at NNLO 
is neither clearly preferred nor rejected by the data shown in Fig. [STJ with some bins favouring 
NLO and some NNLO. Future more precise data may be more discriminating. 

As outlined in Section [U the heavy flavour treatment acquires some model dependence at 
NNLO, in the detailed form of the coefficient functions, which are approximated using 

threshold and small-x limits. However, for charm production this term only appears with all 
scales fixed at = ml in the fit to data (only depending on for < m^), and is a 
small constant correction, which becomes increasingly less important at high Q^. The model 
dependence is mainly due to the treatment of the small-x limit. Reasonable variations of the 
parameters controlling this leads to variations in F2'^{x, Q^) of approximately ±0.03 for x = 10~^ 
or ±0.06 for x = 10~^, and rather less than this for higher x, the variation becoming smaller 
than ±0.01 at x = 0.001, and quickly dying away above this. This is enough to move the fit 
quality for F2'^{x, Q^) data close to that at NLO in one direction, or to move the fit quality 
in the = 2 GeV^ bin close to the best possible in the other direction. Because this is a 
correction to F2{x,Q'^) which is constant in Q^, this uncertainty has very little impact on the 
gluon distribution. The only real uncertainty it leads to is in the input for the sea distribution, 
which would change to compensate the Q^-constant change in F2{x, Q"^). The possible change is 
up to 4-5% for X < 10~^, becoming much less at higher x. This is a lot smaller than the input 
uncertainty shown in Fig.[T9l(a). The relative uncertainty from this source would quickly become 
much smaller at high Q^, and would be a negligible influence on the uncertainty at = 100 
GeV^ shown in Fig. [191 

In Fig. [HHlwe show the predictions for the beauty structure function F|''(x, Q^) compared to 



the published HI data [43|, [4j]. Clearly at all orders the comparison is good. It is also clear that 
these data currently provide no useful constraint on the gluon distribution. Since much of the 
data is in the vicinity of = ml, or not much higher, more precise data will be sensitive to the 
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Figure 37: The charm structure function, F^'^lxyQ'^), compared to data from HI 4l|-|4J] and 
ZEUS 45l-l47l|. The bins used correspond to the data points in Ref. 46|; the other data points 
have been shifted to these values using the NNLO predictions. 
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Figure 38: Predictions for the beauty structure function, compared to data from 

HI HQ. 



91 



details of threshold contributions and in particular to the value of nif, rather than to the gluon. 



9.2 Constraints on intrinsic charm 

As explained in Section HI our variable flavour number scheme allows for the possibility of 
including intrinsic charm with no theoretical problems. However, we simply choose not to 
include it in the default fit since we have no reason to believe it is anything other than a very 
small component of the total charm distribution, and none of the data in our global fit provide 
any constraint. The nonperturbatively generated contribution to the charm distribution is in 
general of 0{AQQ^/ml), but may be enhanced at large x [i^, whereas all data on F2'^{x,Q'^) 
from HERA are at small x. 

In this section we use the only data which currently do show any sensitivity, the EMC data on 



F2'^{x, Q"^) |l65l |. In Fig. [39] we show the comparison of the EMC charm data to our predictions 
at NLO and NNLO. In general the comparison to data is quite poor, tending to overshoot at the 
lower values of x and Q^. This is worse at NNLO than at NLO because the NNLO threshold 
corrections cause a significant increase. Hence, at fixed NLO and NNLO with rric = 1.4 GeV 
there is some room for intrinsic charm at x > 0.2, but much more tendency to overshoot data at 
smaller x. This latter problem must be addressed before drawing conclusions about the intrinsic 
charm contribution. When considering this comparison it is pertinent to remember that when 
calculating heavy flavour structure functions we use W"^ = 4m^ as the threshold, whereas in 
reality we need to produce mesons, i.e. is really required to be a little greater. In practice 
this is a higher- twist effect, but the EMC data are not so far from threshold, and are much 
more sensitive to these details than HERA data. Hence, we consider a comparison where we 
replace by m^(l + A^/m^) in the threshold dependent parts of coefficient functions, where 
is a binding energy which we take to be A = 0.2 GeV. We make no change in the off-threshold 
parts of coefficient functions, nor in PDF definitions. The result is labelled as NLO MT in 
Fig. [391 It clearly matches lower EMC data better, and the prediction is slightly smaller 
at higher x (even with this modification NNLO predictions would generally be too large). The 
curve labelled NLO MT + IC includes a 0.3% integrated number density contribution of intrinsic 



charm (and the same for anticharm) using the parameterisation in Ref. [98| with the upper limit 
on X modified to the correct threshold value. We note that this is only ^ 1/10 of the upper limit 
obtained by CTEQ j99|. This seems to be about the maximum that can be included even with 
nonperturbative suppression of standard fixed-order contributions, and even then only at NLO. 
A much larger contribution would overshoot most data at x > 0.10 Hence, if the EMC data 
are to be believed, there is no room for a very sizeable intrinsic charm contribution, although 
there is a suggestion of some deficit at the highest x values. 



^^Note that 0.3% was presented as the most hkely value in 1983 using very old PDFs and only the LO 
perturbative contribution ^166] . 
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Figure 39: The comparison of the EMC charm data [l65| to our predictions at NLO and NNLO. 
MT stands for the modified threshold approach and IC stands for inclusion of intrinsic charm. 
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10 Description of low-energy Drell— Yan data 



As in recent an alyse s |l8l . |2l| we fit to the E866/NuSea Drell- Yan dilepton production data 
in pp collisions 107l |. now corrected for electromagnetic radiative corrections; moreover we do 
not include the pd collision data as we observe significant systematic differences between our 
predictions and these results. The Drell- Yan process is known up to NNLO for the rapidity 



(y) distribution [167[. However, the data are binned in the Feynman-x variable, xp = Xi — X2- 
Previous MRST fits used LO kinematics {pt = 0) when transforming between the two variables. 
However, the mean pt of each data bin is available, and hence we now use the exact kinematics: 



y 



\n(H + VJpTi) , where H = — ^ . (87) 



At LO a good fit to the data is impossible while simultaneously fitting to structure function data 
which probe the same range of x and hard scale. This is because there is a large difference in the 
size of the corrections to the respective coefficient functions when going from the spacelike to 
the timelike regime, i.e. there is a factor of 1 + {as{M)/7[)CpTT'^ /2 for Drell- Yan production at 
NLO. Since this is missing at LO, the normalisations of the two calculated quantities cannot be 
made consistent with the structure function data and Drell- Yan data at the same time. Previous 
MRST fits used a constant factor of 1.3 to enhance the LO Drell- Yan cross section, but we have 
improved this and now use a i^T-factor consisting of the 7r^-enhanced term. At NLO and NNLO 
we isolate the as dependence from the i^'-factor by writing: 



ir^O(M) = 1 



«s(M)\ CpTi 



.2 



vr / 2 



K^^^{M,y) = l+[^^] CiM,y), (89) 



TT 



ir^^LO(M,y) = l+(^^) D{M,y)+(^^^y E{M,y). (90) 
We calculate the C-, D- and E-factors with the dyrap program [167l | using NLO or NNLO PDFs 



at the respective order. This method, first used in Ref. [21|], is a considerable improvement on 
previous MRST fits where only the overall K{M,y) was calculated for a fixed value of as- 
(Additionally, NNLO fits prior to Ref. used only NLO iC-factors.) 

The quality of the fit to the E866/NuSea data is illustrated in Fig. |40]at NNLO. As seen in 
Table[2l the fit quality is very similar at LO and NLO and the figures of data against theory would 
be similar to that at NNLO. As one can see, there is no systematic failure in the comparison 
between data and theory — the relatively large value is due to a large scatter in the data 
points. The only exception is that in the lowest mass bin the data lie somewhat above the theory 
at low rapidities, although low M is where additional theoretical corrections are most likely to 
play some role and also where there may be complications in the mass reconstructi on f rom the 



proximity to the J/ip and T resonances. Additionally, unlike alternative studies 168| . we do 



not see any significant inconsistency between the E866/NuSea pp data and the rest of the global 
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Ratio of E866/NuSea pp Drell-Yan data (x 1.087) to 
MSTW 2008 NNLO PDF fit, = 237 for 184 pts. 
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Figure 40: The ratio of E866/NuSea pp Drell-Yan data points to the predictions of the NNLO 
PDF fit as a function of rapidity for different mass bins. 
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fit. Our for tlie remaining data only decreases by a few units when these Drell-Yan data are 
removed, i.e. by much less than the tolerance on any of our eigenvectors. We do note, however, 
that the positive corrections of order 10% to the cross section at NNLO require the data to move 
up to a normalisation of 1.087, i.e. slightly outside the quoted uncertaint y of 6.5%. 

As in previous analyses we include the E866/NuSea pd/pp ratio data 108l |. which constrains 
the d — u difference. The experimental extraction of the pd/pp ratio by E866/NuSea used 
a larger data set than the individual cross section measurements, although there was some 
overlap. Moreover, the analysis was of a different form using different cuts and averages over 
much larger data subsets, and therefore had little sensitivity to those lower statistics regions 
where the discrepancies with the pd data occur, in particular at the highest xp values. We 
have improved the theoretical treatment of the pd/pp ratio compared to previous analyses by 
calculating the separate i^'-factors for pp and pd collisions. This has an effect of at most 2%, 
and does aid the quality of the global fit very slightly. There is no significant change in our 
results for the difference d — u, though the uncertainty is reduced a little compared to previous 
analyses, as seen in Fig. [181 
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11 W and Z boson production at the Tevatron 



W and Z boson production at high-energy hadron colhders has traditionally provided important 
constraints on parton distributions. While low-energy Drell-Yan production directly probes the 
quark sea, W and Z production (via qq V, with V = W^, Z) at proton-antiproton colliders 
are mainly sensitive to the valence u and d quarks, in different combinations to structure function 
measurements in deep-inelastic scattering. 

Because of the overall experimental luminosity uncertainty (of order ±6% at the Tevatron), 
the total W and Z cross section measurements do not provide competitive constraints on the 
PDFs, although they do provide an important cross-check on the parton distributions and on 
the hard-scattering framework (see Section [15] below). The shape of the W and Z rapidity 
distributions, however, is sensitive to the shape of the parton distributions, and does provide a 
powerful constraint. We have for many years used CDF data on the W ^ iu charge asymme- 



try [51[ in our global fits, and in the following subsection we study the impact of recent CDF and 
D0 measurements in the updated MSTW fit. There are also nowadays very precise Tevatron 
data on the Z rapidity distribution, which provides complementary information to the W —>■ iu 
charge asymmetry measurement. We include these data in the global fit for the first time, and 
study their consequences. 

As f or th e Drell-Yan process, the W and Z rapidity distributions have been calculated to 



NNLO |l67l |. and so we are able to perform consistent fits at LO, NLO and NNLO. 



11.1 Description of Tevatron W asymmetry data 

The W charge asymmetry at the Tevatron is defined by 

_ da{W^)/dyw - da{W~)/dyw u{xi)d{x2) - d{xi)u{x2) , , 

w[yw) - da(H^+)/d|/H^ + daiW-)/dyw uixi)d{x2) + ^(^^^(xs) ' ^ ' 

where Xi^2 = / y/s) exp{±y]y) . This description in terms of the PDFs is valid at leading 
order when sea-quark contributions are neglected. In practice, it is usually the lepton charge 
asymmetry which is measured, defined in a similar way as 

^_ dain/dve-daii-)/dv, 
da{n/dv^ + da{i-)/dv,' 

where r/i is the pseudorapidity of the charged lepton. Defining the emission angle of the charged 
lepton relative to the proton beam in the W rest frame by cos^ 6* = 1 — 4E^/M^ leads to 

1 / 1 + cos ^* \ 
y, = y^ + -\n{- . (93) 



2 Vl-cos^\ 

Care must be taken when using a valence-quark-only approximation to the lepton asymme- 
try, since near the edges of phase space, i.e. cos^^* ~ ±1, sea-quark contributions can become 
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important. In fact, neglecting overall factors, the numerator of fl92p can be approximated bj@ 

u{xi)d{x2){l — cos6*Y + d{xi)u{x2){l + cos9*Y — d{xi)u{x2){l + cos9*y — u{xi)d{x2){l — cos9*Y, 

(94) 

where for cos ^* f» 1 the leading {du) sea-sea contribution for W'^ production is enhanced relative 
to the valence-valence contribution by the large (1 + cos^*)^ term arising from the V — A decay 
to leptons. Hence if the lepton charge asymmetry is measured in different bins of lepton Et, 
the sensitivity to antiquark PDFs is greatest at lower Et- For example, if Et = 25 GeV then 
cosO* = ±0.78 and, taking the positive value, then from the antiquark contribution will 
be enhanced relative to the quark contribution by a factor (1 + cos^*)^/(l — cos 6'*)^ = 68. 
In this case, in the forward direction, leptons from produced from quarks will be mainly 
from bosons with rapidity yw = y£ + 1.05 (taking cos 6'* = —0.78), while leptons from 
produced from antiquarks will be mainly from bosons with rapidity yw = Z/^ — 1-05 (taking 
cosO* = +0.78). At high 7]£ ^ y£ (not too close to the kinematic limit) this can lead to the latter 
being a significant proportion, or even dominating. Hence, in addition to measuring the overall 
shapes of the quark distributions, the lepton charge asymmetry also probes the separation into 
valence and sea quarks. In practice, this is particularly so for the less well-constrained down 
quark. 

In the global fit, we of course include all quark and gluon contributions to the appropriate 
order in perturbation theory. We calculate ii"-factors for the differential cross sections da{i^) / drj^ 
using the FEWZ code [tsI with LO, NLO or NNLO PDFs at the respective order@ We also 
include the effect of a finite W width in these i^-factorsQ 

In Fig.miwe show the effect of a finite W width, NLO corrections and NNLL py'-resummation 
on the lepton pseudorapidity distribution when using a common set of parton distributions 
(MRST 2004 NLO). Then in Fig. |12] we show the corresponding predictions for the lepton 
charge asymmetry. In the lower E^ bin, finite W width and higher-order effects are important 
only in the forward direction, meaning that these features have a significant effect on the asym- 
metry shown in Fig.HU (Recall that CP invariance requires that da{i~^)/dr]i = da{i~)/d{—rii).) 
In the upper E^ bin, these features are almost symmetric about r/g = 0, meaning that the effect 
on the asymmetry is small. From Fig. |12] we see that at NLO the effect on the asymmetry of 
the finite W width is smaller than the experimental uncertainties, but is not insignificant, and 
hence we include this. The further effect of resummations is generally smaller still, so given that 
our fit is meant to correspond to PDFs at a fixed perturbative order we do not include the effect 
of resummations in our i^-factors. 



^^Contributions involving s, c and b quarks are not displayed, as these are very small at Tevatron energies. 

^■^We find that the VEGAS Monte Carlo integration in fewz does not converge at NNLO for lepton production 
for bin widths of Arji = 0.2 to sufficient accuracy of 0(1%) wit hin a rea sonable run time of C'(months). This 
problem has been noticed by other groups (see, for example, Refs. 169l . 17d |). Since the effect of NNLO corrections 



is expected to cancel almost entirely when calculating the asymmetry Aijig), we instead use NLO isT-factors 
calculated with NNLO PDFs. We have confirmed that the NNLO corrections to the W asymmetry (as a function 



of the W rapidity) are extremely small; see Fig. 11 of Ref. [167 |. 

^''The future use of APPLGRID ^171] may eliminate the need to resort to calculating LO W and Z cross 
sections multiplied by fixed X-factors during PDF fits, although the latter method should still be a very good 
approximat ion . 
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Figure 41: da/drif,+ i^-factors for CDF kinematic cuts 48| for various methods of calculation 
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Figure 42: W ^ eu charge asymmetry from CDF [48| for various methods of calculation 
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CDF data on lepton charge asymmetry from W^ev decays 
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Figure 43: A comparison of theory and data for W ^ ev charge asymmetry from CDF |4 
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D0 data on lepton charge asymmetry from W^^iv decays 
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Figure 44: A comparison of theory and data for W ^ fdv charge asymmetry from D0 [49 



Figure |13] shows the fit to the most recent CDF W ev lepton charge asymmetry data [48 . 
The data are presented in two bins of lepton Et- In each case the data are reasonably well 
fitted at both NLO and NNLO. Also shown is the NLO fit with the antiquark contributions 
removed. The difference is much more significant in the lower Et bin, as expected. In the higher 
Et bin there appears to be some evidence for a systematic discrepancy between the fit and 
the data, particularly at high |?7e| (these data seem to be described slightly better by CTEQ6.6 
PDFs). This could in principle be removed by increasing the d quark distribution in the x ~ 0.1 
region, thereby enhancing the negative d{xi) contribution in the numerator in flMl) . but this 
would produce a tension with DIS structure function data. The overall impact on the ci-quark 
distribution will be discussed further below. 

The corresponding fits to the D0 W 
shown in Fig 



49 



are 



lepton charge asymmetry data from Ref. 
The data are presented in a single lepton > 20 GeV bin, and the antiquark 
contributions are therefore sizeable. Some systematic discrepancy at high \ri^\ is again apparent. 
Note that there are more recent D0 electron data S^l than the muon data shown in Fig. HH 
However when we include these D0 electron data in the global analysis we are unable to obtain a 
good quality fit, with significant tension between the new D0 electron asymmetry data and the 
DIS structure function data (both proton and deuterium) and with low-mass Drell-Yan data. 
Furthermore, we notice considerable tensio n be tween the D0 electron asymmetry data ||50|] and 
the preliminary CDF W asymmetry data 172],^ Pending further investigation of these issues, 
we have decided not to include either the more recent D0 data 50| or the CDF W asymmetry 



25 



The CDF W asymmetry data have since been pubhshed [17 
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data 172. 173 in the current fit. 



11.2 Description of Tevatron Z rapidity distribution data 

The integrated luminosity for Run II at the Tevatron has allowed a high-precision measurement 
of the Z rapidity distribution. Analogously to the low-energy Drell-Yan distribution d^a/dMdy 
discussed in Section [TOl this provides a constraint on the (mainly quark) PDFs over a broad 
range of a; at ~ M|. The mix of quark distributions probed is approximately 

0.29m(xi)m(x2) + 0.37c;(xi)rf(x2), (95) 

in contrast to the (4mm + d(£) combination probed in Drell-Yan production. The Z rapidity 
distribution is therefore particularly sensitive to the d distribution, providing complementary 
information to the W charge asymmetry discussed above. 

The theoretical calculations are performed using the same form of i^-factors as for the low- 
mass Drell-Yan data fl88l) . fl89]) and fl90|) . with the NLO and NNLO coefficients calculated exactly 
using the VRAP program [l67 ]. 

In Fig. US] we show the results of the NLO and NNLO fits compared to the Zj'^* rapidity 
shape distribution measured in the mass range 71 < Mee < 111 GeV by D0 using 0.4 fb~^ of Run 
II data jiil . In Fig. |16] we show the results of the NNLO fit compared to the preliminary Zj^* 
rapidity distribution measured in the mass range 66 < Mee < 116 GeV by CDF using 2.1 fb~^ 



of Run II data [52] • The data are well fitted, particularly when the systematic uncertainties 
are taken into account. There is evidence of a slight excess of data over theory at high rapidity 
{y > 2.2 ^ xi > 0.4 and X2 < 0.005). 

During the fit the Z rapidity distribution was calculated at LO in the narrow- width approx- 
imation (Mee = Mz) neglecting the 7* contribution, then applying the appropriate i^-factors 
to account for higher-order QCD corrections. To be more precise, the invariant mass range 
66 < Mee < 116 GeV should be explicitly integrated over in the theory prediction allowing for 
a finite Z width and including the 7* contribution. This more precise treatment has the effect 
of lowering the theoretical prediction by less than 1% over most of the rapidity range, rising 
above unity at \y\ ~ 2, then rising more rapidly as the kinematic limit is approached at still 
larger rapidities. Refitting the PDFs with this additional factor included gives an increase in 
the fitted CDF normalisation by 1% and a better description of the CDF Z/7* rapidity data 
at large rapidity, with smaller shifts needed in the correlated systematic uncertainties. For the 
NLO fit, the foi' the 29 CDF data points improves from 49 to 42, while the fo^' the 28 
D0 data points improves from 19 to 17. However, the refitted PDFs are almost unchanged, 
with the differences very small compared to the quoted uncertainties. Since the CDF data are 
still unpublished and may change before publication, we will update our fits in the near future 
using the final Z/'-f* data and accounting for the precise Mee range and 7* contribution in the 
theoretical calculation. 
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Z/y* rapidity sliape distribution from D0 
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Figure 45: The comparison of data and theory for the Z/'-f* rapidity shape distribution from 

D0 ii. 
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Z/y* rapidity distribution from CDF 

MSTW 2008 NNLO PDF fit, yj- = 50 for 29 points 
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Figure 46: The comparison of data and theory for the Zj^* rapidity distribution from CDF 52 
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Figure 47: The impact of new data on the down valence quark distribution at NNLO. 

11.3 Impact on the down valence quark distribution 

The W and Z data at the Tevatron are mainly sensitive to the up and down quark distributions 
of the proton. However, the up quark distribution is already well constrained by proton structure 
function data where it appears charge-weighted. Therefore, in practice, the Tevatron W and Z 
data mainly constrain the down quark distribution. In Fig. H?] we compare the down valence 
distribution at NNLO with that from the MRST 2006 NNLO set. Note that we now make a 
better choice of parameters in the input dy distribution ([7]) to allow for more flexibility for the 
eigenvector PDF sets, which generally gives larger uncertainties than previously, despite the 
presence of more constraining data. Clearly there is a significant change in the shape, though 
the new and previous versions generally agree within the 90% CL. limits. This change is mainly 
due to the inclusion of new Tevatron data, though there is also some influence from the new 
neutrino structure function data. The most significant change is the increase in d^lxjQ"^) for 
X ~ 0.3, seen more clearly in Fig. [61T b). which is influenced by the lepton asymmetry data. 
From the number sum rule this has to be compensated for by a decrease elsewhere, which is 
mainly for x ~ 0.07. 



106 



12 Description of jet data 



Knowledge of the high-x gluon distribution is important for calculations of both signal and 
background in new physics searches at the LHC The main direct constraint comes from Tevatron 
data on inclusive jet production. This is well illustrated by looking at the analyses of the Run I 
Tevatron jet data. 



12.1 Description of Run I Tevatron jet data 



Historically, early CDF Run I inclusive jet data at central rapidity jl74j indicated an excess in 
the high-ii^T jet spectrum using NLO QCD with the available parton distributions of the time. 
This discrepancy was initially taken as a possible indication of quark compositeness. It was found 
to be "impossible" ^ to describe simultaneously both the CDF jet data and the DIS structure 
function data with a common set of parton distributions having a physically motivated input 
form. However, it was found that the data could be accommodated in a glob al fit with an input 



gluon parameterisation with a somewhat contrived shoulder at high x |l75l |. Subsequentl y, D0 



Run I jet data became available, which measured the Et distribution at a range of rapidities |176 



NLO global analyses were able to fit all the jet data with physically reasonable forms, though 
there was some tension between the description of the high-i?^ CDF and D0 data. This "saga" 
indicated the need to provide PDFs with associated uncertainties so that it can easily be seen 
in which regions of x the PDFs are rel ative l y co nstrained or unconstrained. (These issues are 
further discussed in, for example, Refs. |l77l . Il78[ |.) 

The Tevat ron Run I jet data have been routinely included in the default global analyses 
since CTEQ4 [iT^ and MRST 2001 
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The MRST 2004 analysis [18| showed that a physical 
parameterisation of the input gluon distribution could be obtained by taking the input form 
(|T0|) in the DIS scheme then transforming to the conventional MS scheme. The resultant MS 
parameterisation was found to automatically give a shoulder-l ike form at large x, which produced 
a better description of the Tevatron Run I inclusive jet data |176l . Il79| than the previous MRST 
2001 PDFs. 



At NLO, the cross section for inclusive jet production at the Tevatron can be written as 



app = alifil) [a^^ + as{fil)a^^°] ® fa/p{xa, ® fb/p{xb, 



(96) 



=q,9 



Usually, multidimensional (Monte Carlo) phase space integration is used to evaluate fl96l) with 
cuts on, for example, the transverse energy Et and pseudorapidity t] of the produced jets. 
However, this can take hours or days of CPU time, so it is impractical to include jet data in 
a PDF fit in this way. The old (approximate) solution used by MRST was to calculate "i^'- 
factors", defined as CTpp'^/app, for a given set of PDFs. Then six pseudogluon "data" points at 
fijp = 2000 GeV^ (and a value of Aqcd corresponding to asi^M^) ~ 0.115) were inferred that 
would correspond to the best app needed to describe the data. These six pseudogluon points 
(with a free normalisation) and the Aqcd point were included in the MRST fit rather than the 
original jet data points. The comparison to the actual jet data points, including the inclusion 
of the correlated systematic uncertainties, was only made after the fit, though it was checked 
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that the changes induced in the pseudodata were quantitatively similar to those in the final 
comparison to the jet data. The CTEQ group do include the jet data points explicitly in their 
fit, but the jet cross section is still calculated at LO and a NLO /^-factor applied. 

A much improved solution, which we adopt here, is to interpolate the PDFs and as so that 
they can be factorised from aab in ( l96l) . Then the convolution integrals of ( l96l) are replaced with 
multiplication. The phase space integration then only needs to be done once for a given set of 
kinematic cuts, with the result store d in a grid for later use during a PDF fit. This technique 



has pr eviou sly been used by HI |l8l| to determine the gluon density from DIS jet data, and by 
ZEUS 140l | in a PDF fit including jet data from DIS and photoproduction. Grids corresponding 
to the kinematic cuts used in published Tevatron and HERA iet data are provided by the 

(Another project with 



FASTNLO project [75 



76, 77 



which uses the nlojetH — \- calculations 
similar aims, which has a lso b een applied to calculate W and Z production using mcfm 
is the APPLGRID project jl7l| .) 



12.2 Description of Run II Tevatron jet data 



Recently, Tevatron Run II inclusive jet data from both CDF SJ, l55|] and D0 [56j] have become 
available with much higher statistics and much better control over the corrections necessary 
to go from the "raw" data to measurements at the parton "jet" level. For example, D0 now 
claim a relative uncertainty on the jet pt of less than 2%. Our present global analysis uses 
all the available Run II jet data, including the complete information on correlated systematic 
uncertainties in the definition fl38l) . Due to the increased reliability, and statistics, of the Run 
II jet data, we now omit the Run I data from the global fit. 

Another change is that we now use a renormalisation and factorisation scale Hr = fJ^p = Pt, 
rather than fiR = fxp = pt/2@ The latter scale appears to be very convergent for central 
rapidities, but becomes unstable for larger y. For example, the NLO contribution is some 60- 
70% of that at LO for y ~ 2.5. In comparison, the scale choice fiR = fir = Pt leads to NLO 
corrections of ~ 20-30% for all y (except at the highest y and pt, where they are smaller) and 
hence, while not as convergent near y = 0, the scale choice ^r = fip = pt is reasonably stable 
for all rapidities of the data@ 

In Fig.SHlwe show the comparison of the NLO fit with CDF Run II inclusive jet data using the 
kx jet algorithm ^]. There is also an analysis of the same data using the cone-based Midpoint 



algorithm [55|, but the kr algorithm data [5J] were published and implemented in our analysis 
first, therefore only the kr algorithm data are included in the final global fit. Furthermore, there 
is a theoretical preference for using the k^ jet algorithm rather than a seeded cone algorithm, 
since the latter is not formally infrared safe, which in practice means that the perturbative 
calculation is only accurate to a given order of the expansion (NLO for the Midpoint algorithm). 
In Fig. |49]we show the comparison of the NLO fit to D0 Run II inclusive jet data using a cone 



jet algorithm 56||. We achieve an excellent description of both data sets shown in Figs. HHl and 



^^Note that the Run II data are binned in terms of transverse momentum px and rapidity y, while the Run I 
data were binned in terms of transverse energy Et and pseudorapidity rj, although Et = Pt and rj ^ y for the 
massless jets used in Run I, but not for the massive jets used in Run II. 

^^We thank M. Wobisch for bringing this to our attention. 
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CDF Run II inclusive jet data, ^ = 56 for 76 pts 
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Figure 48: Data/theory ratio for the CDF Run II data obtained using the fcr jet clustering 
algorithm Q and the MSTW 2008 NLO fit. 
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D0 Run II inclusive jet data (cone, R = 0.7) 

MSTW 2008 NLO PDF fit (|i = |i = p^^^^), = 114 for 110 pts. 
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Figure 49: Data/theory ratio for the D0 Run II data [56( and the MSTW 2008 NLO fit 
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Note that the data points are allowed to shift by the correlated systematic errors according 
to the definition fl38l) . Recall from Table [3] that the CDF Run II inclusive jet data are forced 
to share the same normalisation as the CDF Run II Z rapidity distribution [iil] due to the 
common luminosity uncertainty. In practice, the shared normalisation of 0.96 at NLO (or 0.99 
at NNLO) is determined by the much more constraining Z rapidity distribution. 

In Fig. [50] we compare the predictions using the MSTW 2008 NLO PDFs to the CDF Run 



II data obtained using the Midpoint cone algorithm 55|, which were not included in the fito 



The agreement is generally good, but with the second and third point for 0.7 < \y 



JET I 



and the first 
for 72 points. 



)oint for 1.1 < |?/"^^^| < 1.6, each contributing about 11-13 to the total = 108 
□ These three points appear simply to be outliers. Replacing the CDF Run II 
data measured using the /c^ jet algorithm by these cone algorithm data and refitting improves 
the by 6 and the Xgiobai by 4. The x^ of the D0 Run II cone algorithm data is unchanged 
on refitting. We conclude that there is no significant tension between the three Tevatron Run 
II data sets on inclusive let production |54|-|56| obtained using either the cone [55|, l56| or kr 
jet clustering algorithms 
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12.3 Impact on the gluon distribution 

The impact of the Tevatron Run II jet data on the high-x behaviour of the gluon distribution is 
interesting. These data prefer a smaller high-x gluon distribution than the Run I data. Indeed, 
we find that it is no longer necessary to take the input gluon parameterisation fIlOp in the DIS 
scheme and then transform to the MS scheme, as was done in the MRST 2004 [18] and MRST 



2006 [21| analyses in order to give an enhanced gluon at high x to better fit the Tevatron Run 
I data. Instead, the input gluon form ffTOj) is now taken directly in the MS scheme. 

In Fig. [51] we show the high- a; gluon distribution from the MSTW 2008 NLO analysis com- 
pared to that obtained from a fit with the Tevatron Run II jet data replaced by Tevatron Run 
I jet data, and from a fit without any Tevatron jet data. The fit to Tevatron Run I data uses 



FASTNLO 



75| and includes hadronisation corrections that were previously neglected in all global 
fits, but can reach the level of a (9(10%) correction in the lower Et bins. We manage to fit the 
Tevatron Run I jet data well using the input gluon parameterisation (fTOj) in the MS scheme. 



without needing to perform the DIS— > MS transformation used in the MRST 2004 18| and 



mg 

2i 



MRST 2006 [21| analyses. The 68% C.L. uncertainty bands from the fits to Run I and Run 
II data clearly disagree; the 90% C.L. uncertainty bands also disagree for x > 0.4 (not shown 
here) . 



^^We include the complete 24 (excluding luminosity) sources of correlated systematic uncertainty ^55)]. Note 
that the data tables originally published in Ref. [H^ have received minor corrections, and the corrected data are 
shown in Fig. [50] 

^^We note that the rapidity region 0.7 < |y''^"^| < 1.1 corresponds to a crack between the central and plug 
calorimeters [55l |. 

^°If, instead of tying the normalisation to the CDF Run II Z rapidity distribution [55|, we allow a separate 
normalisation for the CDF Run II inclusive jet data, the preferred normalisation at NLO for the Midpoint 
algorithm data is 1.06, with an improvement in to 91 for 72 points. By contrast, the kx algorithm data prefer 
a normalisation of 0.98 with almost no improvement in compared to the default normalisation of 0.96. 
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CDF Run II inclusive jet data, = 108 for 72 pts 



8 1.8^ 



a 
Q 



1.6 
1.4 
1.2F 
1 

0.8P 

0.6 

0.4 



0.0 < ly'^^^l < 0.1 



.5 



10^ 



(GeV) 



S 1.8 

t 1.6 
n 

n 1.4 

a 

1.2 



1 — 
0.8 r 
0.6 r 
0.4 L 



0.7< Iv'^^^l < 1.1 



o o O o 



« o « o * J c > a a ^ 



10^ 



pJET (GeV) 



o 



1.6 
1.4 
1.2 
1 

0.8 
0.6 
0.4 



L 1.6< |y''^V2.1 

L < 




'-_ oooooooOc^Sl 

. . . 1 





10^ 



(GeV) 



1 1.8 r 

t 1.6 
n 

n 1.4 
O 

1.2- 



1 — 

0.8 r 

0.6 r 
0.4 L 



0.1 < Iv'^^^l < 0.7 



, * * ^ ' o * c? J ' ■ 



o o § 



o o o o 



10^ 



(GeV) 



S 1.8 

t 1.6 
n 

ra 1.4 
Q 

1.2 



0.6 
0.4 



L 1.1 < ly-^^^l < 1.6 

. ..... n.o'fft 






. . .1 





10^ 



pJET (GeV) 



Midpoint: R = 0.7, f =0.75 

merge 

MSTW 2008 NLO PDF fit 

^ Without systematic 
uncertainties 



Withi systematic 
uncertainties 



Figure 50: Data/theory ratio for the CDF Run II data obtained using the cone-based Midpoint 
jet clustering algorithm jssi and the MSTW 2008 NLO fit. These data are not included in the 
fit. 



112 



(a) 

Gluon distribution at Q^ = W'^ GeV^ 




Figure 51: The impact of the Tevatron jet data on the high-x gluon distribution at NLO. 
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UDri 179 
(33 pts.) 


D{ol 176 
(90 pts.) 


LDbLL[kT) |54J 
(76 pts.) 


D(Z)11 |56J 
(110 pts.) 


A 2 

A.non-jet 

(2513 pts.) 


as{Mi) 


53 
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64 


117 





0.1197 


51 


48 


132 


180 


9 


0.1214 


56 


110 


56 


114 


2 


0.1202 


53 


85 


68 


117 


1 


0.1204 



Table 8: values for Tevatron data on inclusive jet production from various global fits. The 
large bold values indicate the data sets explicitly included in each fit. The "MSTW 2008" fit 
includes only the Run II data. We also indicate the increase in the for data sets other 
than Tevatron inclusive jets when these data are added to the fit. Finally, the value of a5(M|) 
is given for each fit. All fits are carried out at NLO with a scale choice of fiR = fip = px- 



We also compare in Fig. [5T] to the high-x gluon distribution of the ZEUS 2005 Jets analy- 
sis 



140|, which is even smaller than our fit without the Tevatron jet data. This is likely to at 
least partially be due to the omission of fixed-target data in the ZEUS 2005 Jets analysis. In 
a global analysis, where as{M'^) takes a value of 0.117-0.120, the fixed-target data constrain 
the high-x gluon distribution to be moderately large, though in general the constraint is corre- 
lated with the value of as{M'^). The ZEUS 2005 Jets analysis also has a more re strict ive input 
gluon parameterisation than our distribution. Note that although the ZEUS fit 140| included 
additional data on dijets in photoproduction, which are not included in our fit, these data were 
found only to reduce the uncertainties without changing the central values, so do not provide 
an explanation for the discrepancy. 

In Table IH] we show the values for the Tevatron inclusive jet production data from various 
global fits. The four lines in the table are fits (i) without Tevatron jet data, (ii) with only Run 

I data, (iii) with only Run II data, and (iv) with both Run I and Run II data included. The 
CDF Run I data 179l | exist only for a central rapidity bin, and hence are less sensitive to the 
gluon distribution than the other data sets in Table [Hi The P^r data point when fitting only 
Run I data is significantly larger than 1 for CDF, suggesting unusually large fluctuations, and 
significantly smaller than 1 for D0, suggesting an overestimation of errors. The per data 
point when fitting only Run II data is much closer to 1 for both CDF and D0, suggesting that 
the Run II data are more reliable than the Run I data. The tension between the Run I and Run 

II data is evident from Table [S], i.e. the fit to only Run I data gives a poor description of Run II 
data, and vice versa. The global for the fit without Tevatron jet data is 2371 for 2513 points. 
When including the Run I data it is 9 worse for the same data points, but when including the 
Run II data it is only 2 worse. The as values are also slightly different: the fit without Tevatron 
jets has as{M'^) = 0.1197 compared to the fit with Run I data which has a5(M|) = 0.1214 and 
the standard fit with Run II data which has a5(M|) = 0.1202. The Run II jet data are therefore 
slightly more consistent with the other data sets included in the global fit than the Run I jet 
data, and lead to a significantly lower value of as, which influences the gluon distribution. 

In Fig. [52r a) we show that the inclusion of Tevatron jet data does reduce the fractional 
uncertainty on the high-x gluon, albeit not dramatically. This rather disappointing result is 
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Figure 52: (a) The impact of the Tevatron jet data on the fractional uncertainty of the high-x 
gluon distribution at NLO. (b) The impact on the high-x gluon distribution of replacing CDF 
Run II jet data obtained with the kT jet algorithm 5J] by data obtained with the Midpoint cone 
algorithm 55|]. Also shown is the effect of taking an alternative scale choice. 
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due to the large number of sizeable correlated systematic uncertainties on the jet data, which 
are free to adjust to best fit the data according to the definition fl38|) . There is therefore a 
considerable trade-off between these systematic shifts and the input gluon parameters. 

The high-x gluon distribution from the fit using CDF Run II data obtained with the Midpoint 



cone algorithm [55j instead of the kx algorithm [5^ lies slightly below the central value of the 
MSTW 2008 NLO fit, but is well within the uncertainty band, lying slightly above the gluon from 
the fit without Tevatron jet data; see Fig. [52](b). (The value of ^^(Ml) = 0.1196.) This fact 
confirms the statement made in Section [12.21 that the two CDF Run II data sets are consistent. 
In Fig. [52r b) we also show the effect of taking a scale choice of fiR = fip = Pt/'2 rather than our 
default choice of fiR = fip = Pt for Tevatron jets. The lower scale choice results in a smaller 
high-x gluon distribution, which is within the 90% C.L. uncertainty band (but outside the 68% 
C.L. uncertainty band, not shown here). The value of as{M^) = 0.1192 is also smaller with 
the lower scale choice. The quality of the description to either of the CDF Run II data sets is 
similar with either scale choice, but the of D0 Run II data is ~10 worse with the lower 
scale choice. Therefore, overall, the Tevatron Run II jet data prefer our default scale choice of 
IJ'R = f^F = Pt- By contrast, when fitting only to Tevatron Run I jet data, the of the D0 
Run I data is ~10 better with = = PtI'^i i-e. the default scale choice used in previous 
MRST and CTEQ analyses. Finally, in Fig. [521(b) we also show the high-x gluon distribution 
obtained from a simultaneous fit to Run I and Run II jet data, i.e. the last row in Table [H 
Unsurprisingly, it lies in-between the gluon distributions obtained from the two separate fits to 
either Run I or Run II jet data alone. Note, however, that all gluon distributions shown in 
Fig. [521(b) are obtained from fits with slightly different values of as-, and hence are not really 
directly comparable. 

A comparison with earlier gluon distributions is shown in Fig. l53l The plots show that the 
inclusion of the Run I Tevatron jet data in the MRST 2004 [l8j and CTEQ6.6 [87] analyses 
led to a significantly harder gluon than the present analysis (MSTW 2008) which is based on 
the Run II Tevatron jet data. (For the MRST 2004 set, where eigenvector PDF sets are not 
available, we assume that the fractional PDF uncertainties are the same as for the MRST 2001 
set.) In fact, with the Run II data, we are almost back to the "soft" MRST 2001 gluon, which 
did not describe the Tevatron Run I jet data nearly as well as the MRST 2004 analysis. 

As explained in Section [STTl we include the Tevatron Run II jet data in the NNLO analysis on 
the basis that the NNLO theoretical correction is very likely to be small and a smooth function 
of Pt, and to be less than the size of the correlated systematic uncertainties on the data. As can 
be implied from Figs. [51] and [52], the inclusion of these data does not influence the central value 
of the high- a; gluon distribution at all significantly, but does help constrain the uncertainty. 

12.4 Jet production at HERA 



In Figs. [5l] and [55] we show the comparison of the NLO fit to ZEUS data |57l . 1580 on inclusive jet 
production, and in Fig.l56]we show the comparison of the NLO fit to HI data [5 91]. Jet data have 
been included previously by the ZEUS Collaboration in their own PDF analysis ll^. HI have 
similarly included inclusive jet data in their determination of diffractive PDFs |l82| . However, 
we find that the HERA jet data have little influence on either the central values or the PDF 
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Figure 53: Comparison of the current high-x gluon distribution at NLO with that from the 
previous MRST and CTEQ [87| analyses. 
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Figure 54: Data/theory ratio for the ZEUS 96-97 inclusive jet data [57[ and the MSTW 2008 
NLO fit. 
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Figure 55: Data/theory ratio for the ZEUS 98-00 inclusive jet data [58( and the MSTW 2008 
NLO fit. 
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Figure 56: Data/theory ratio for the HI 99-00 inclusive jet data [59|] and the MSTW 2008 NLO 
fit. 
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uncertainties: the constraint they provide with the current accuracy on the data being mimicked 
by other data sets in a fully global fit. We do not include the dijet photoproduction data as this 
involves further systematic uncertainties due to the choice of photon PDFs. 

We omit the HERA jet data from the NNLO fit since it is quite probable that the unknown 
NNLO corrections are fairly large. However, the description of these data is excellent if we use 
NLO hard cross sections with NNLO parton distributions, i.e. comparable to the fit we obtain 
at NLO. If the NNLO correction is equivalent to a fairly smooth (and not too large) JC-factor, 
the good fit could likely be maintained by the interplay of data relative to theory using the 
correlated systematic uncertainties. 
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13 Low-x gluon distribution and Fi 



In this section we outline our results for the small-x gluon distribution and the implications 
of a measurement of Fl{x,Q'^). This is, in principle, a very important quantity to compare to 
predictions since it gives an independent test of the gluon distribution at low x to accompany the 
indirect determination from dF2{x,Q'^) / dhiQ"^ , which beyond LO results in a perhaps surpris- 
ingly small, or even negative, gluon distribution. However, it is also a direct test of the success 
of extensions of our theory beyond fixed-order calculations. Until recently we have been limited 
to consistency checks on the relationship between F2{x, Q^) and Fl{x, Q"^) at high y, where both 
contribute to the total cross section a"(x, Q"^) = F2{x, Q^)—y'^/{l + {l—yy)FL{x, Q'^). Extracting 
Fl data from this requires an extrapolation in y making some theory assumptions. It is more 
useful to fit directly to a{x, Q"^) where there is a turnover at the highest y. This was discussed 
in Section I7.H where we saw that although the full turnover in Q^) implied by the data 
was not seen at NLO, or even in full for the central fit at NNLO, the quality of fit to the data 
was acceptable given the large uncertainties. Clearly the precision of such studies is limited, 
and can be affected by systematics, e.g. the photoproduction background uncertainty. However, 
using the final low-energy run at HERA, a direct measu rement o f Fl{x, Q^) at HERA has been 



possible, and some results have already been presented |183l . Il84l . 

The prediction for Fl{x,Q'^) is mainly determined by the form of the gluon distribution 
extracted from the fit. This is shown in Fig. [571 There is poor stability as one changes the 
perturbative order of the calculation, particularly at small x and Q^. At small x there is a 
large order- by-order change in the splitting functions, particularly Pgg, and this leads to a large 
variation in the gluon. At NLO a small-x divergence appears for the first time, making the 
evolution of the quarks much quicker, and hence the NLO gluon is dramatically smaller at small 
X. Most extractions find a NLO gluon which is small and perhaps valence-like at < 2 GeV^ 



and very small x. However, the gluon parameterisation we have used since Ref. allows the 
possibility for the gluon distribution to become negative, and for the best fit it chooses to do so 
for X < 10~^ at = 2 GeV^. It is frequently claimed that higher twists, in particular gluon 
recombination effects, will be a cure for this, but sim ply a pplying these corrections, on top of 



NLO leading-twist QCD, results in little change [17|, Il85l |. At NNLO there is an additional 
ln(l/a;)/a; divergence, and at the smallest x the NNLO gluon becomes smaller to compensate, 
or in practice more negative for very small values of x. However, it is important to note that 
the gluon beyond leading order does not have a simple probabilistic interpretation. This is 
particularly the case at high orders at small x, where perturbative corrections to cross sections 
can be large, and even more so in unphysically motivated schemes such as the MS factorisation 
scheme. 

The fact that the gluon is not directly physical is well-illustrated by the NNLO 0{al) lon- 



gitudinal coefficient function C^^^g{x) [64, 165|. This has a large positive contribution at small x, 
and this counters the decrease in the small-a: gluon. The predictions for Fi{x,Q'^) at LO, NLO 
and NNLO are shown in Fig. [581 The Fi{x, Q"^) prediction is more stable in going from NLO to 
NNLO than the gluon distribution at small x. Despite the behaviour of the gluon distribution, 
the central value of F^ is positive at NNLO for all x > 10~^ at = 2 GeV^, and is greater 
than at NLO at the smallest x. The uncertainty becomes enormous as x decreases below 10""^, 
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Figure 57: The gluon distribution at LO, NLO and NNLO including the one-sigma PDF uncer- 
tainty bands. 
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Figure 58: Predictions for the longitudinal proton structure function, Fl^x^Q"^), at LO, NLO 
and NNLO includi ng t he one-sigma PDF uncertainty ban ds. Also shown are resummed NLL 



BFKL predictions |144| | and dipole model predictions (186 
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but at the lowest the NLO and NNLO predictions are discrepant in some regions. The LO 
prediction is far larger than either, reflecting the huge correction in the small-x gluon going to 
NLO. In order to more clearly separate the relative effect of different gluon distributions, in 
Fig. [59]we show the predictions for Fl with LO, NLO and NNLO coefficient functions but using 
common NNLO PDFs. There is a positive correction at high x at each order. NLO then results 
in a negative correction at lower x, whereas NNLO gives an additional negative correction at 
intermediate x before becoming positive again at sufficiently small x. However, the value of x at 
which these transitions occur is very sensitive to . As the gluon becomes steeper the regions 
of the gluon distribution probed in the convolution become more local to the value of x, and 
small-x divergences in the coefficient functions become less important, along with the relative 
importance of higher orders decreasing as the coupling becomes weaker. Hence, the ultimate 
positive effect of the ln(l/a;)/x term in the NNLO coefficient function is only seen to increase 
the NNLO result above LO for x > 10~^ in the two lowest bins. 

As the relative lack of stability in fixed-order predictions implies, there are various potentially 
large corrections beyond fixed-order perturbation theory. It is pos sible that there is a large 
higher- twist contributi on fr om renormalons in the quark sector 187l |. and the implications are 



discussed in Refs. |28l . Il88| . Perhaps more significantly, since the small-x NNLO correction is 
itself rather large, even higher orders might be important. There are leading ln(l/x) terms of 
the form: 



XPgg 



4- 



xPqg{x^ 



n-2 



(1/x) and xCi^gi 



X] 



4- 



[l/x). (97) 



A fit which performs a double resummation of leading and next-to-leading ln(l/x) and run- 
nin g cou pling contributions leads to a better fit to small-x data than a conventional perturbative 



fit 144| . The gluon distribution from this resummed fit, defined in a more physical scheme, is 



larger at small x and than NLO or NNLO, and indeed is always positive for 1 GeV^. 



189, 190 



all lead 
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This is reflected also in the prediction for Fl{x,Q'^). Similar approaches 
to rather comparable results for the calculated splitting functions, but only in Ref. 
detailed phenomenological studies taken place. A partially overlapping set of additional correc- 
tions are considered in the dipole picture. As with small-x resummations, this approach can be 
cast in the language of f{x,k^), the unintegrated gluon distribution, which is directly related 
to the dipole-proton cross section. The structure functions are obtained by convoluting this 
dipole cross section with the wave functions for the photon to fluctuate into a quark-antiquark 
pair. Hence, this picture includes some of the resummation effects, but also higher-twist contri- 
butions, and is designed to approach = smoothly. However, it misses quark and higher-x 
contributions. In this framework, higher-twis t correcti ons are not small in either or Fj- sep- 
arately, but largely cancel in F2 = Fl + Ft |l43l . Il9l| . Overall the Fl{x,Q^) predicted in the 



dipole model approach is steeper at small x than fixed-order predictions, and is automatically 
stable at lowest Q^. The general features are rather insensitive to whethe r sat uration effects are 
included in the d ipole cross section. Resummed NLL BFKL predictions 144| and dipole model 
predictions jl86| are additionally shown in Fig. ESI illustrating these features^ 



The first published direct measurement of Fl{x, Q"^) can be found in Ref. 183l | for 12 < < 



90 GeV^, based on data taken in the last few months of HERA running, when the proton beam 
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Figure 59: Predictions for the longitudinal proton structure function, Fii^x^Q"^), using NNLO 
PDFs and LO, NLO and NNLO coefficient functions. 
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corn- 



Figure 60: The longitudinal proton structure function, Fl{x,Q'^), measured by HI 
pared to predictions at LO, NLO and NNLO incl udin g the one-sigma PDF uncerta inty bands. 
Also shown are resummed NLL BFKL predictions 144 1 and dipo le model predictions 18d |. (The 



dependence on the form of the dipole cross section is mild dij.) 



energy was lowered from the nominal value of 920 GeV to values of 460 GeV and 575 GeV. 
In Fig. [60] we show these data compared to our predictions calculated at LO, NLO and NNLO 
using the appropriate parton distributions at each order. The uncertainty bands obtained from 
the 40 alternative eigenvector PDF sets are shown at each order. In the region of the existing 
HI data the NLO and NNLO predictions agree well with the data, while the LO prediction 
overshoots the data at the lower values. We also show in Fig. [60] the predictions at lower 
< 10 GeV^. We extrapolate the values of x accordingly, by fitting the x values of the HI 
data as a power law function of Q"^, obtaining x = (1.09 x 10~^)((5^/Qq)^-^^ with Ql = I GeV^. 
The NLO and NNLO predictions are very similar along this line, but as seen in Fig. [58] this is 
partially accidental. Note that the NNLO uncertainty band at low x and is smaller than 
the NLO uncertainty band due to the lesser sensitivity to small-a; PDFs in the convolution 
integrals at NNLO compared to NLO, because the more divergent coefficient f uncti on samples 
the PDFs further fr om x . Again we also show resummed NLL BFKL predictions jl44| and dipole 
model predictions \lS(^ , and both agree well with the published data. For E [5, 10] GeV^ 
the uncertainty in NLO/NNLO predictions for Fl{x, Q"^) due to the gluon uncertainty increases 
to more than 20%. A good measurement of Fl{x,Q'^) here would automatically improve the 
gluon determination. Resummations and dipole models suggest a higher low-Q^ Fl{x,Q'^) by 
an absolute value of up to 0.15 — well outside the fixed-order uncertainties. So, moreover, a 
good measurement of Fl[x,Q'^) in this region would start to discriminate between theories. 
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14 Comparison with other PDF sets 



In this section we make a comparison with various other PDF sets obtained from analyses of 
a variety of data sets by different groups. Some of the most important features pertaining to 
previous MRST sets have aheady been mentioned, but here we will present a general overview 
of the similarities to, and differences from, other sets, and suggest reasons to explain these 
differences. 

We begin with a historical perspective. In Fig. [6T] we show the current NNLO PDFs com- 



pared to those from MRST 2006 [21[ . The comparison would in general terms be rather similar 
if comparing our current NLO set to the (unpublished) 2006 NLO set, but because the change in 
a5(M|) has been slightly larger at NNLO some differences in PDFs are a little more apparent in 
the NNLO comparison. Comparing with older NLO sets would mean a change in the definition 
of the GM-VFNS, and the consequences of this have already been illustrated in Fig. [31 Beginning 
from the top of Fig. EH we see that the shape of both valence distributions has changed sig- 
nificantly, though there is still compatibility with MRST 2006 within the 90% CL. uncertainty 
bands. This is partially due to the smaller value of ^^(Ml) — the coefficient functions give a 
positive effect at very high x and negative for x ~ 0.1, and this is smaller in 2008 than 2006, so 
the quarks compensate. However, the details, and the large changes in are due to the new 
Tevatron Run II W and Z data and the NuTeV and CHORUS neutrino DIS data. For x < 0.01 
the form of both valence quark distributions is determined mainly by the requirement to satisfy 
the number sum rules. The gluon distribution has changed significantly: the smallness at high-x 
being mainly due to the Tevatron Run II jet data, as discussed in detail in Section [121 This 
is part of the reason for the decrease in 05 (M|), and the weaker coupling and the momentum 
sum rule both contribute to the significantly larger gluon distribution at lower x in the 2008 
set. This feeds through, via evolution, to the sea quarks being higher at small x, though well 
within previous uncertainties. For the gluon and light quarks the uncertainty is generally a little 
larger than for the MRST 2006 set. This is due largely to the inclusion of data normalisation 
uncertainties into the PDF error analysis, but also due to an improvement in the flexibility of 
the parameterisation for the gluon and down valence distribution, and for light sea quarks the 
independent determination of the strange distribution also leads to an increase. The fact that 
we fit the strange distribution directly results in the very large increase in the uncertainty for 
X > 0.001. But as we assume that the input strange has the same shape as x — >• as the light 
quark distributions, the uncertainty becomes the same as for u and d at very small x and high 



In Fig. [Mlwe show the NLO PDFs at = 10^ GeV^ compared to those from CTEQ6.6 [87 



The CTEQ6.6 analysis uses a very similar total set of data, but is not as up-to-date, in par- 
ticular not including any Tevatron Run II data or the neutrino DIS data sets from NuTeV and 
CHORUS. Overall, there is good agreement between the two sets, with very few regions where 
the 90% C.L. limits do not overlap. It is striking that the difference between the MSTW 2008 
and CTEQ6.6 valence quarks is rather similar to that between the MSTW 2008 and MRST 2006 
distributions. This suggests that the difference is largely due to the new NuTeV, CHORUS and 
Tevatron Run II data included in our analysis. The strange quark distribution is rather different, 
though the 90% C.L. uncertainty bands overlap. This is due to the different assumptions made 
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Figure 61: MSTW 2008 NNLO PDFs compared to MRST 2006 NNLO PDFs. 
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Up valence distribution at = 10" GeV^ 
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Down valence distribution at = 10" GeV^ 
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Up quark distribution at Q = 10 GeV^ 
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Figure 62: MSTW 2008 NLO PDFs compared to CTEQ6.6 NLO PDFs. 
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in parameterising the strange distribution, and was addressed in Section I6.5[ As in previous 
comparisons the gluons have a systematic difference. Since MSTW use the ffexible parameteri- 
sation which allows the distribution to be negative at very small x at input (and at 1.69 GeV^, 
the CTEQ input scale) while the CTEQ input is valencelike, the latter is bigger at the smallest 
X values. The gluon in the CTEQ6.6 set is determined at high x by Run I Tevatron jet data, 
so is higher in this region, and from the momentum sum rule must be smaller at intermediate 
x ~ 0.1 to compensate. However, the greater similarity in the heavy flavour treatments now 
employed perhaps leads to better agreement in the intermediate x range than previously. The 
comparison of the u and d distributions for x < 0.1, i.e. outside the valence region, very much 
reflects the gluon difference. Finally, we note that the MSTW 2008 uncertainty bands are not as 
small compared to CTEQ bands as in previous comparisons, despite the addition of more data 
and on average a smaller tolerance for the former. This is a reflection of the data normalisation 
uncertainties and the improved flexibility in parameterisations. 



In Fig. [63] we show the NLO PDFs compared to those from Alekhin 1391 ]. As one can 
see there are many significant differences. The Alekhin NLO PDFs are obtained from a fit to 
structure function data only. An update i nclud ing some Drell-Yan data, and in particular the 
E866/NuSea pd/pp data, was made in Ref. 168| . but only NNLO sets were made available. Since 
NLO sets are more widely used, and there are some systematic differences in the approaches at 
NNLO, we instead compare to the NLO set. There is clearly a significant discrepancy in quark 
distributions everywhere, despite the fact that largely the same neutral-current DIS data are 
fit. This highlights the importance of fitting neutrino DIS data and Drell-Yan and vector boson 
hadroproduction data to obtain the quark flavour decomposition. The absence of Tevatron jet 
data plays a part in the gluon distribution in the Alekhin set being much smaller at x > 0.3. The 
smaller value of a5(M|) = 0.117 also contributes to differences, as does the fact that the PDFs 
correspond to the ZM-VFNS. However, these certainly do not explain the major systematic 
differences. 

Finally we compare to the very recent NNPDFl.O set [l37| . which uses proton and deuterium 
target DIS data, as well as CHORUS data (without nuclear corrections), again in the ZM-VFNS. 
This approach relaxes the restriction in uncertainty due to fixed parameterisation limitations, 
though in practice relies on a very large number of parameters. It also proceeds by making a 
Monte Carlo sample of the distribution of experimental data by generating a large number of 
replicas of data centred on each data point with full inclusion of the information from errors 
and their correlations. Each replica is used to generate a PDF set and the mean for a quantity 
obtained from averaging, and uncertainties from standard deviations. The data are split into 
training and validation sets for each replica, and the fo^^ oiis is monitored while the other 
minimised, thus avoiding over-complicating the input PDFs by stopping when the fit to the 
validation sets stop improving. As with the Alekhin set there are some significant systematic 
differences between these and the MSTW 2008 sets. However, the use of some neutrino DIS data 
and the fact that a fixed asiM^) = 0.119 is used in NNPDFl.O, together with a more flexible 
input parameterisation compared to Alekhin, means that there are less dramatic variations. 
Some of the most significant variations are shown in Fig. [M] where we compare our up valence 
and up antiquark distribution to those from NNPDFl.O [137]. Some significant differences are 
also observed in the total up quark distribution and the small-x gluon at high Q^. From Fig. [M] 
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(a) 



Up valence distribution at = 10" GeV^ 



MSTW 2008 NLO (68% C.L.) 

Alekhin 2002 NLO 




(b) 



Down valence distribution at = 10" GeV^ 




(d) 



Up quarl< distribution at 0=10" GeV^ 




Down quarl< distribution at =10" GeV^ 

MSTW 2008 NLO (68% C.L.) 




(f) 



Strange quark distribution at = 10" GeV^ 




Gluon distribution at = 10" GeV^ 




Figure 63: MSTW 2008 NLO PDFs compared to Alekhin NLO PDFs. 
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Figure 64: MSTW 2008 NLO (a) up valence quark and (b) high-x u distribution compared to 
those from NNPDFl.O. 



we see there are sizeable differences in the central values of the up valence distribution, even 
where the sea is very small, perhaps due to the fact that relative data set normalisations are fitted 
in our approach but not in the NNPDFl.O analysis. The difference in the valence distribution is 
also undoubtedly due to the lack of data distinguishing between valence and sea quarks. Indeed, 
the most striking disagreement is in the up antiquark distribution at very large x. In the global 
fit this would be constrained to be positive by Drell-Yan data. Some differences in central 
values is not a great surprise but it is more surprising that the error propagation used in the 
NNPDFl.O analysis does not completely encompass the MSTW 2008 values. Perhaps this is a 
reflection of the fact that some data sets used in the fit are incompatible, and discrepancies will 
sometimes occur, as we found comparing our benchmark fit and standard fit for the high-x gluon 
in Section [6l6l Since it is the uncertainties that are currently most interesting for the NNPDFl.O 
set we concentrate our comparison on this. In Fig. [65] we compare our fractional uncertainties to 
those from NNPDFl.O 137|. For total quark and gluon distributions the NNPDF uncertainties 



are only a little larger, and it must be remembered that we have many more data sets in the fit. 
For the strange quark the fact that we fit directly, whereas in NNPDFl.O the strange distribution 
at input is a fixed fraction of the total light sea, actually leads to our uncertainty being larger 
for 0.01 < a; < 0.2. The increased flexibility in the NNPDF parameterisation does lead to a 
larger uncertainty at high x and small x for valence distributions, though we can argue that 
we use more constraining data for 0.01 < x < 0.7, and have less flexibility from the sum rule 
outside this range. The NNPDFl.O valence quark uncertainty narrows in the region x ~ 0.025 
at = 10^ GeV^ (and at x ~ 0.06 at = 2 GeV^) which we do not understand, as this is 
just below t he lo wer limit of the BCDMS data, which is one of the best constraints. An update, 
NNPDFl.l 159( 1 . has recently been reported. The uncertainties are larger in this analysis due 
both to a randomisation of the preprocessing exponents and due to extra freedom introduced in 
the strangeness sector. Since there are no dimuon data used to constrain the strange PDF, and 
very little constraint from any other data set, the uncertainty on s{x, Q'^) + s{x, Q^) in NNPDFl.l 
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a) 



Up valence distribution at = 10" GeV^ 




(b) 



Down valence distribution at = 10" GeV^ 




(d) 



Up quarl< distribution at = 10" 

MSTW 2008 NLO (68% C.L.) 



GeV^ 




Down quarl< distribution at Q = 10" GeV^ 




Strange quark distribution at = 10" GeV^ 




(f) 



Gluon distribution at = 10" 




Figure 65: MSTW 2008 NLO PDF uncertainties compared to those from NNPDFl.O. 
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is huge. This feeds into the other PDFs, with uncertainties increasing by factors of three for 
many quark combinations, and it is difficult to tell what effect the randomisation alone has. This 
illustrates how important the direct constraint on the strange PDF, and some theoretical input in 
the parameterisation, can be when obtaining PDFs with genuinely representative uncertainties. 
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15 W and Z total cross section predictions 



In this section we present and discuss predictions for the total cross sections for W and Z 
production multiphed by leptonic branching ratios at the Tevatron and LHC. To begin with we 
consider W = + W~ production, and subsequently consider the separate and W~ total 
cross sections at the LHC. (At the Tevatron, the and W~ total cross sections are equal.) 

To avoid confusion, we first define exactly what we mean by the total W and Z cross sections. 
We assume a narrow-width approximation, in which the W and Z are treated as on-shell parti- 
cles, of mass Mw and Mz respectively, in the calculation of the production cross sections. These 
cross sections are then multiplied by the appropriate leptonic branching ratios. We assume that 
the cross sections measured at the leptonic level over a restricted region of phase space can be 
corrected back to this "theoretical" quantity. We present predictions at LO , NL O and NNLO 



in perturbative QCD, using a program based on th at of van Neerven et al. [193| incorporating 



subsequent corrections from Harlander and Kilgore [19J|. We use LO electroweak perturbation 
theory, with the qqW and qqZ couplings defined by 

gl, = GfM^/V2, gl = GfMIV2. (98) 

The electroweak mixing angle sin^ 9^ determines the relevant strengths of the vector and axial- 
vector Z couphngs to the quarks. Since our cross section predictions depend to some extent 
on the elect roweak parameters, we use the latest (2008) values from the Particle Data Group 



(PDG) [i32j- 

We sum over u, d, s, c, b quarks in the initial state, using the 2008 PDG (unitarity 



constrained) values of the CKM weak mixing matrix elements \ Vqqi \ in the case of production. 

The MSTW 2008 predictions for the total W and Z cross sections, including the one-sigma 
PDF uncertainties, are given for the Tevatron and LHC in Table |9l It is seen that the NNLO 
PDF uncertainties are around 2% at both the Tevatron and LHC. The uncertainty coming from 
neglected higher-order QCD corrections beyond NNLO, estimated by allowing the factorisation 
and renormalisation scales to vary so that the ratios uf/M w.z, fJ^R/Mw,z and iif/ij^r all lie in 



the range between 0.5 and 2 (as in, for example, Ref. IQSf ). is around ±0.5% at the Tevatron 
and ±1% at the LHC. Uncertainties not well accounted for by scale variations, in particular the 
effect of small-x resummations, variations in GM-VFNS and changes in as{M'^) may all lead to 
further uncertainties. These are difficult to quantify at NNLO, but at the LHC may be 1-2% 
in each case. For the Tevatron only the last, i.e. as'(Ml), is likely to be significant. In Fig. [66 
we compare the MSTW 2008 NLO and NNL O pr edictions for the W and Z total cross sections 



to Tevatron data from CDF |l96l | and D0 [197|. The agreement is good. The experimental 



error, which is dominated by the ~ ±6% uncertainty on the machine luminosity measurement, 
is significantly larger than the uncertainty on the theoretical predictions. 

For both theory and experiment, the uncertainty in the W and Z cross section ratio is much 
smaller than on the individual cross sections. This is illustrated in Figs. [67] and [68], which 
show two-dimensional plots of the MSTW 2008 NLO and NNLO W cross sections against the 
corresponding Z cross sections, including the error ellipses calculated using f[M|) and f[55|) with 
the one-sigma error setslll We compare to the central values of cross sections calculated with 



^^As noted in Ref. [87|, the error ellipse calculated with e.g. one-sigma PDFs corresponds to a confidence level 
smaller than one-sigma in two dimensions 131 . 132| . 
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Table 9: Predictions for W = + W and Z total cross sections at the Tevatron and LHC, 
including the one-sigma PDF uncertainties, and their ratio Rwz- We take fJ,R = fip = Mw,z- 
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Figure 66: W = + W~ and Z total cross sections compared to Tevatron data |l96l . Il97l |. 
where the dashed lines show the one-sigma PDF uncertainties on the NNLO predictions. The 
luminosity uncertainty of ~ ±6% is not included on the experimental data points. 



137 



W and Z total cross sections at the Tevatron 
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Figure 67: W = + W and Z total cross sections at the Tevatron {y/s = 1.96 TeV). The 
error eUipses, calculated using the one-sigma error sets, are shown for the MSTW 2008 NLO and 
NNLO PDFs. We show the CDF Run II measurement of the ratio Rwz = 10.84 ± 0.15(stat.) ± 
0.14(syst.) jl96 | with statistical and systematic uncertainties added in quadrature to get the 
overall uncertainty band. 
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W and Z total cross sections at the LHC 
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Figure 68: W = + W and Z total cross sections at the LHC {^/s = 14 TeV). The error 
eUipses, calculated using the one-sigma error sets, are shown for the MSTW 2008 NLO and 
NNLO PDFs. 
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Ratio to MSTW 2008 
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Table 10: Ratios of predictions for W = + and Z total cross sections at the Tevatron 
{y/s = 1.96 TeV) and LHC {y/s = 14 TeV) calculated using the central values of different PDF 
sets, with respect to those from the MSTW 2008 sets. 



a variety of other PDF sets. For the LHC predictions in Fig. [68], we show lines of constant 
Rwz = Bii, ■ aw / Bii- az, and on the corresponding Tevatron plot in Fig . [671 we overlay a recent 
measurement of Rwz = 10.84 ± 0.15(stat.) ± 0.14(syst.) from CDF [igj- 



The ratios of predictions for W and Z cross sections calculated using the central values of 
different PDF sets, with respect to those from the MSTW 2008 sets, are given in Table [TOl The 
difference between the predictions of the various PDF sets can be understood by considering the 
values of the corresponding PDFs in the relevant x and region. We have already noted that 
the relevant PDFs in the MRST 2006 and MSTW 2008 sets are very similar, see Fig. [6ll and 
this is reflected in the similarity of the W and Z cross sections. We have also noted that a more 
careful treatment of the heavy quark (c and b) distributions at NNLO in the MRST 2006 and 
MSTW 2008 sets gave rise to a significant increase in the light quark distributions at high Q^, 
and this explains the increase in the W and Z cross sections compared to MRST 2004o The 
CTEQ6.6 NLO predictions are close to those of MSTW 2008 NLO at the Tevatron, but some 
2% higher at the LHC. This can be traced back to the more rapid evolution of the CTEQ quarks 
at small x, driven by a (larger) gluon that, unlike MSTW 2008, is constrained to be positive at 
Qq, see Fig. [621 The Alekhin NLO and NNLO predictions for the ratio Rwz at the LHC are 
larger than the predictions from all other PDF sets, due to the small strange quark distribution 
in the Alekhin PDF sets; see Section. [TH 

We have so far considered total W = + W~ production, but of course at the LHC 
(Tw+ 7^ cind the ratio R± = aw+f^w- is also of interest from the PDF perspective. 

Assuming a diagonal CKM matrix we have, at leading order, 

^ _ u{xi)d{x2) + c(xi)g(x2) + (1^2) 
^ ^ d{xi)u{x2) + s{xi)c{x2) + {I ^ 2)' 

For central production, which dominates the total cross sections, we have Xi ~ X2 ~ Mw / a/s and 
Q ~ Mw- In the limit ^/s — >■ cxo, R± —>■ 1 since the quark and antiquark distributions become 
equal in this limit. The fact that R± 7^ 1 at LHC is a reflection of the residual differences. 



similar effect was observed at NLO in the evolution of CTEQ6.1 to CTEQ6.5 and CTEQ6.6 [8J,[87| 
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Table 11: Predictions for and W total cross sections at the LHC, including the one-sigma 
PDF uncertainties, and their ratio i?±. We take (j^r = fip = Mw- 

particularly between u and it and between d and d, due to the valence quarks. Fig. [69] shows 
various predictions for the and W~ total cross sections at the LHC. Lines of constant R± 
are also superimposed. The MSTW 2008 predictions are listed in Table [TTl 

Comparing Fig. ED with the corresponding Fig. ESI for W and Z production, we see that 
there is less correlation between and W~ than between W and Z. This is because the 
combinations of u and d quark and antiquark distributions probed in W = + W~ and Z 
production are more similar than in and W~ separately. It is also interesting that the 
prediction for R± has decreased significantly in going from MRST 2006 to MSTW 2008, due to 
a change in the u/d ratio resulting from the addition of the new Tevatron W and Z data and 
neutrino DIS data, as discussed in Section [TTl 

Note that unlike Rwz, which additionally depends on electroweak parameter and branching 
ratio values, the overwhelmingly dominant uncertainties in the theoretical predictions for R± 
are those due to the PDFs. The experimental measurement should also be very precise, since it 
is simply a matter of comparing the number of and i~ events in a sample of W ^ iu events. 
From Fig. [^Hl we see that a measurement of R± with an error of less than 1% at the LHC will 
further constrain the parton distributions, particularly the u/d ratio. 

In Fig. [To] we show the ratio of the parton luminosities 

|^ = ^/^/a(x,Mi)/,(r/x,Ml), r=^ (100) 

at y/s = 10 TeV compared to y/s = 14 TeV at the LHC for ab = gg,'Yliq=u bl'ly using the 
MSTW 2008 NLO PDFs including the one-sigma uncertainty bands. The PDF uncertainty on 
the ratio of the W total cross section at 10 TeV compared to 14 TeV is of order ±0.3%, and 
similarly for the ratio of Z total cross sections. The PDF uncertainty is so small because it 
is only sensitive to the slope (rather than absolute values) of the PDFs in the region between 
X = 0.006 and x = 0.01 (at central rapidity). Other theoretical uncertainties, such as the choice 
of electroweak parameters, should also cancel out in this ratio. Experimentally, the accuracy of 
the measurement of this ratio will depend on whether the relative machine luminosity at 10 TeV 
and 14 TeV can be measured with high precision. 
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and W total cross sections at the LHC 
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Figure 69: and W total cross sections at the LHC {^/s = 14 TeV). The error elhpses, 
calculated using the one-sigma error sets, are shown for the MSTW 2008 NLO and NNLO 
PDFs. 
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Figure 70: Ratio of parton luminosities at i/i = 10 TeV compared to y/s = 14 TeV at the LHC, 
using the MSTW 2008 NLO PDFs including the one-sigma uncertainty bands. 
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16 Conclusions 



We have presented the new MSTW 2008 parton distribution functions (PDFs) at leading-order 
(LO), next-to-leading order (NLO) and next-to-next-leading order (NNLO). In each case we 
include eigenvector sets spanning the uncertainty of the PDFs defined at a fixed value of a5(M|). 
These sets are a major update to all our previously available sets, and incorporate the maximum 
amount of information from deep-inelastic scattering (DIS) and other hard-scattering data. We 
have also introduced a number of new theoretical refinements to produce what are arguably 
the most complete and reliable sets available. The PDFs are specifically designed for use in 
precision cross section predictions and uncertainties at the LHC The grids and interpolation 
code (in Fortran, C++ and Mathematica formats) are available from hepforge [13], and are 



also made available as part of the lhapdf package |198l ] from version 5.7.0. 



The principal new data sets included in the analysis are: 

• neutrino structure function {F2 and xF^) data from NuTeV 37|] and CHORUS 



38 



neutrino dimuon cross sections from CCFR and NuTeV 39|, which constrain the strange 



quarks and antiquarks, allowing these to be fit directly for the first time, instead of using 
the previous input assumption that s = s = n{u + d) /2 with k ^ 0.5; 



the Z rapidity distribution [52, |53|] and the lepton charge asymmetry from W decays [48 



49| measured at the Tevatron Run II, both of which constrain mainly the down quark 



distribution; 



the inclusive iet production data from the Tevatron Run II [54, 156|, which prefer a smaller 
high-x gluoio than the previous Run I data; 



• updated HERA data on |4l|-|47l| and on inclusive jet production in DIS |57H59[ . 

In addition to including new data, there are a number of important developments in the 
fitting procedure. 

• Our treatment of heavy flavours, carried out in a general-mass variable flavour number 
scheme (GM-VFNS), is the most complete available, particularly at NNLO where the 



improved treatment, first implemented in the MRST 2006 analysis [2l|], has a significant 



quantitative effect. At NLO we estimate the possible theoretical uncertainty on the parton 
distributions due to variations in the specific choice of GM-VFNS to be around 2%, though 
this should decrease at NNLO. 

We improve and extend our input parameterisation, not only introducing two new free 
parameters for each of s + s and s — s, but using one more in the determination of the 
uncertainty on the gluon distribution, and making a better choice for the down valence 
quark distribution. 



^^The smaller gluon distribution at high x has important consequences. For example, the predicted Tevatron 
cross section for production via gluon-gluon fusion of a Standard Model Higgs boson with mass Mh = 170 GeV 



decreases by almost 15% compared to the resuh using the previous MRST 2006 NNLO PDFs [199l 1200 1 (for 
implications, see Ref. j20li ]). although this is also partly due to the smaller value of as for MSTW 2008 NNLO. 
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We have designed and implemented a new treatment of PDF uncertainties, doing away with 
the conventional notion of an overall fixed tolerance Axgiohni — ^0 (^^^ ^ 90% confidence- 
level limit) and replacing it with a new "dynamical tolerance" in which the tolerance for 
each independent eigenvector direction is determined from the condition that all data sets 
should be described within their 90% (or 68%) confidence-level limits. In the same context, 
we have also included uncertainties from data set normalisations for the first time. This 
can have an important effect on the PDF uncertainties, particularly for light quarks in the 
range of x and where they are needed for collider phenomenology. It also means that, 
despite the presence of more constraining data in the fit, the PDF uncertainties have got 
slightly bigger in some cases. 

We have replaced the unconventional choice of as used in the MRST fits, which is not 
strictly correct at NNLO, with the exact solution of the renormalisation group equation 
(El), as used in recent pubhc DGLAP evolution codes. This has enabled our evolution code 
to be checked against the PEGASUS 68] and hoppet [69'] codes for the first time. 



We use the fastnlo package [75| for NLO (and for Tevatron jets also for approximate 
NNLO) calculations of jet cross sections directly in the fit. We have also improved our 
calculation of i^-factors f or lo w-mass Drell-Yan lepton pair and W and Z boson production 



using the dyrap /vrap |l67l ] and FEWZ [73] codes 



We have fitted a wide variety of data (~ 2700 data points in all) and overall the quality of the 
NLO and NNLO fits is similar and is perfectly acceptable, with x^/^pts. ~ 1 for almost all data 
sets fitted. As has been noticed in previous studies, however, the quality of the LO fit is poor, 
and there are a number of examples where the data really do require significant NLO corrections. 
We therefore do not recommend that the LO set be used, for example, in leading-order parton- 



shower Monte Carlo event generators. A modified LO set (see Refs. 66|, l67|) suitable for this 
purpose will be produced in a subsequent study. 

Each of our best-fit PDF sets comes with a corresponding set of 40 error PDFs (compared 
to 30 in previous MRST versions), which can be used to estimate the PDF uncertainty on 
any physical quantity. It is important to remember that these PDF uncertainties come from the 
experimental errors on input data only. They do not include additional theory errors, coming, for 
example, from missing higher-order terms in the perturbation series. This issue will be addressed 
in a future study. While the convergence of the series is an important and outstanding issue for 
some quantities studied, for example, small x structure functions, it does not appear to be a 
problem in other cases, for example, (Tw,z at the hadron colliders. 

In subsequent studies we will return to some of the more theoretical issues, including the 
dependence of the fit on the value of and the heavy quark (c^fe) masses, and the impact on the 
fits of including 0{a) QED corrections, as was done in Ref. [l9|]. We also intend to investigate 
W, Z and Drell-Yan lepton pair production at the LHC in more detail, with a particular focus 
on extracting information on the parton distributions at very small x. 

In the meantime, these distributions are obtained from the most up-to-date and complete 
global fit yet performed, and we believe they are the most reliable sets available as we enter the 
LHC era. 
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